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CHAPTER  1:  Introduction 


1.1  Overview 

Detonations  of  energetic  materials  can  lead  to  damage  of  infrastructure  and  loss  of 
life.  Detonations  can  occur  when  the  sensitivity,  of  a  generally  chemically  unstable  energetic 
material,  to  a  thermo-mechanical  stimulus  is  large,  which  would  effectively  reduce  the 
required  energy  to  initiate  a  detonation.  Therefore,  understanding  and  controlling  the 
sensitivity  to  detonation  of  energetic  materials,  through  hot  spot  fonnation,  is  critical. 

Energetic  materials  generally  are  comprised  of  a  crystalline  solid  that  is  encapsulated 
in  a  polymer  binder,  where  the  crystalline  solid  contains  the  expendable  energy  that  is 
released  during  detonation.  Energetic  crystals  commonly  used  are  PETN  (Pentaerythritol- 
tetranitrate),  HMX  (cyclotetramethylene-tetranitramine),  and  RDX 

(Cyclotrimethylenetrinitramine)  [1].  The  combination  of  an  easily  defonnable  polymer  and  a 
generally  unstable  crystalline  material  can  provide  desirable  thermo-mechanical  behavior. 
The  combination  of  polymer  and  crystal,  a  plastic  explosive,  furthermore,  results  in  an 
energetic  material  that  can  be  molded,  machined,  or  otherwise  tailored  for  desired 
applications.  The  polymer  within  an  energetic  material  is  also  needed  due  to  the  lack  of  self¬ 
adhesion  of  the  crystalline  energetic  grains. 

Thermo-mechanical  stimulation,  due  to  defects,  such  as  dislocations,  voids,  and 
cracks  can  cause  energy  to  localize,  which  gives  rise  to  hot  spot  formation  within  energetic 
materials  [2-6].  A  hot  spot,  which  has  achieved  a  critical  size  and  temperature,  causes  the 
energetic  material  to  deflagrate  and  subsequently,  under  specific  conditions,  lead  to 
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detonation  [6].  Many  microstructural  mechanisms  and  characteristics  have  been  suggested  as 
the  underlying  mechanisms  in  understanding  how  hot  spots  form  and  evolve.  Some  of  the 
dominant  observed  microstructural  mechanisms  are  voids  (both  inter-crystalline  and  intra- 
crystalline)  [2,  7,  8],  shear  banding  [9],  and  dislocation  pile-ups  [5,  10],  Studies  by  Borne, 
Mory,  et  al.  [7]  and  Qui,  Stepanov,  et  al.  [11]  also  clearly  indicate  that  as  the  porosity 
volume  fraction  increases,  the  sensitivity  to  shock  initiation  increases. 

The  myriad  heterogeneous  microstructural  aspects  that  span  multiple  physical  scales 
associated  with  energetic  materials,  such  as  grain  boundaries  (GBs),  energetic  crystals, 
polymer  binder  composition,  defects,  and  void  distributions  are  further  complicated  by 
dynamic  pressure  loading  conditions.  The  hot  spot  micromechanics  model  developed  by  Wu 
and  Huang  [12]  attempts  to  predict  the  formation  of  hot  spots  in  drop  weight  experiments. 
This  model  describes  the  microstructure  as  spherical  domains  of  energetic  crystals,  where  the 
grain  morphology  and  defects  were  smeared  out.  Whitworth  and  Maw[13]  developed  a  ID 
Lagrangian  multi-material  model  to  predict  temperature  rises  in  a  single  material  with  one 
void.  Several  molecular  dynamics  simulations  have  also  been  undertaken  to  study  energetic 
materials  [14-19].  However,  defects  such  as  large  voids  and  GBs,  that  span  the  micron  scale 
are  computationally  difficult  to  represent  in  molecular  dynamic  simulations  over  physically 
realistic  time  scales. 

Panchadhara  and  Gonthier  [20]  utilized  a  mesoscale  analysis  to  predict  the  effects  of 
volumetric  and  surface  dissipation  in  granular  HMX.  Their  predictions  indicate  that 
temperature  increases  in  the  material  were  mainly  due  to  plastic  work  heating,  and  the  peak 
temperatures  observed  were  due  to  frictional  heating.  Baer  [21]  developed  a  model  to 
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understand  the  highly  fluctuating  stress  states  in  energetic  materials  through  using  multiple 
stress  jump  states.  A  Lagrangian  finite  element  framework  has  been  developed  that  accounts 
for  polymer  binder  relaxation,  frictional  heating,  and  interfacial  debonding  and  determined 
that  the  fonnation  of  hot  spots  was  due  to  frictional  heating  and  viscous  dissipation  in  the 
binder  [22,  23]. 

Dynamic  behavior  can  lead  to  microstructural  failure  mechanisms,  such  as  void 
coalescence,  crack  nucleation,  and  subsequent  crack  propagation[24-28].  Identifying  and 
understanding  how  cracks  nucleate  and  propagate  is  essential  to  delineating  dissipative 
energy  and  failure  pathways.  In  crystalline-amorphous  systems,  cracks  propagate  in  the 
crystalline  phase  and  typically  tenninate  at  the  crystalline-amorphous  interface,  and  there  is 
generally  a  material  competition  between  dynamic  crack  propagation  and  plasticity.  [29,  30] 

Experimental  work  by  Bouma  et  al.  [24]has  shown  that  energetic  RDX  crystals 
contain  cracks  and  voids  due  to  processing.  Palmer  and  Field  [26]  showed  experimentally 
that  under  compressive  loading  of  HMX  energetic  crystals,  twinning  occurs,  and  this 
typically  precedes  crystal  fracture.  Sharma  et  al.  [27]  has  shown  that  RDX  energetic  crystals 
undergo  large  scale  nano-fracturing  that  breaks  up  crystals  into  nano-crystallites  with  sizes 
ranging  from  500  to  20  nm  for  dynamic  shock  loading  conditions.  Modeling  and 
experimental  work  by  Bouma  et  al.  [25]  has  also  shown,  for  dynamically  defonned  energetic 
materials,  that  fracture  and  crystal  debonding  are  major  failure  mechanisms.  Smith  and 
Thorpe  [28]  showed  that  fracture  in  the  energetic  material,  Composition  B,  occurred  by 
cracking  in  the  RDX  phase  followed  by  cracking  in  the  TNT  matrix  phase.  Experimental 
characterization  work  by  Voevodin  and  Zab inski  [31]  has  also  indicated  that  dynamic 
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fracture  of  crystalline-amorphous  composites  occurs  due  to  cracks  nucleating  in  the 
crystalline  phase,  and  then  tenninating  at  the  crystalline-amorphous  boundaries.  Zhang  et  al. 
[30]  also  demonstrated  that  in  a  crystalline  Cu  and  amorphous  system  that  cracks  nucleate 
and  propagate  in  the  crystalline  phase  and  terminate  at  the  amorphous  phase  boundary.  Wang 
et  al.  [29]  showed  that  fracture  in  a  ductile  crystalline-amorphous  nano-laminate  occurs  by 
cracking  in  the  crystalline  region  where  cracks  terminate  at  the  amorphous  phase  boundary. 

1.2  General  Research  Objectives  and  Approach 

The  objectives  of  this  research  were  to  identify  the  microstructural  defects  and 

mechanisms  that  were  associated  with  hot  spot  formation  and  local  failure  in  energetic 
materials,  and  to  qualify  the  temporal  evolution  of  these  failure  mechanisms.  This  allows  for 
the  identification  of  microstructural  defects  and  mechanisms  such  as  voids,  grain  boundaries, 
grain  morphology,  grain  misorientation,  dislocation  densities,  inelastic  defonnation,  crack 
nucleation,  crack  propagation,  and  interfacial  effects  in  energetic  materials  to  be  associated 
with  hot  spot  formation  and  local  and  global  failure  in  the  material. 

The  effects  of  voids,  high  dynamic  pressure  loading  conditions,  high  dynamic  strain 
rate  loading  conditions,  grain  boundary  misorientations,  crack  nucleation  and  propagation, 
and  viscous  flow  behavior  in  polymer  binder  materials  were  investigated  by  means  of  a 
crystalline  plasticity  finite  element  method  which  had  been  tailored  to  include  the  effects  of 
thennal  decomposition  in  energetic  materials. 

The  microstructural  mechanisms  associated  with  hot  spot  fonnation  specific  to  the 
energetic  material,  RDX-PCTFE,  were  detennined.  However,  the  effects  of  the  PCTFE 
binder  were  studied  using  a  hypo-elastic  constitutive  law,  and  therefore  neglected  the  effects 
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of  viscous  deformation  in  polymer  binders,  which  are  typical  in  energetic  materials. 
Therefore,  the  effects  microstructural  defects  were  investigated  under  the  influence  of  a 
viscous  binder,  estane,  and  the  mechanisms  that  led  to  hot  spot  formation  were  realized. 

Next,  the  effects  of  crack  nucleation  in  an  aggregate  were  studied  to  determine  how  intra- 
crystalline  porosity  affects  crack  nucleation.  The  mechanisms  for  crack  nucleation  on  the 
microstructural  scale  were  obtained  for  local  failure  of  the  material  aggregate.  Next, 
energetic  material  interfaces  were  studied  to  understand  how  interfaces  affect  crack 
propagation  and  the  mechanisms  by  which  crack  arrest  occurs  in  energetic  materials. 

1.3  Dissertation  Organization 

This  dissertation  is  outlined  as  follows:  Chapter  2  presents  the  crystalline  plasticity 
and  dislocation  density  formulation,  the  finite  viscoelasticity  formulation,  and  thenno- 
mechanical  coupling.  Chapter  3  presents  the  numerical  method  used  for  obtaining  quantities 
used  in  the  formulation  presented  in  Chapter  2.  Chapter  4  contains  the  results  from  an 
investigation  into  the  effects  of  microstructural  defects  on  hot  spot  fonnation  in  RDX-PCTFE 
energetic  aggregates.  Chapter  5  covers  the  investigation  into  the  effects  of  a  polymer  binder 
on  hot  spot  formation  mechanisms  for  an  RDX-Estane  aggregate.  Chapter  6  outlines  results 
from  an  investigation  into  the  crack  nucleation  behavior  and  mechanisms  of  RDX-Estane 
energetic  aggregates.  Chapter  7  outlines  results  on  how  interfaces  in  energetic  materials  are 
paramount  to  crack  arrest  and  how  GB  misorientations  affects  crack  propagation. 
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CHAPTER  2:  Crystalline  Plasticity  Formulation,  Finite  Viscoelasticity  Formulation, 
and  Thermo-Mechanical  Coupling 


2.1  Multiple  Slip  Crystal  Plasticity  Kinematics 

The  constitutive  formulation  for  rate-dependent  multiple-slip  crystalline  plasticity, 

coupled  to  evolutionary  equations  for  the  dislocation  densities,  will  be  outlined  below  [32, 
33], 


The  velocity  gradient  tensor,  Ltj,  is  calculated  from  the  defonnation  gradient  as 


L  =f.f:\ 

V  y  V 


(2.1) 


It  is  assumed  that  the  velocity  gradient  can  be  additively  decomposed  into  elastic  and  plastic 


parts,  Le..  and  Lp...  These  will  be  further  decomposed  into  the  symmetric  deformation  rate 

e  p  e  p 

tensors,  Z)  and  Z). ,  and  antisymmetric  spin  tensors,  Wlf  and  Wtj  ,  as 


e  p 

Dy=D,  +  Dy, 


w^w'  +  w*. 


(2.2) 


The  plastic  parts  are  related  to  the  crystallographic  slip  rates  as 


V  9  ? 


j jrP 

WV  =0),  r  , 


(2.3) 


where  a  is  summed  over  all  slip-systems,  and  Py  and  are  the  symmetric  and 


antisymmetric  parts  of  the  Schmid  tensor,  defined  in  terms  of  the  slip  planes  and  directions  as 


PV)  =  !(>)„(«)  +s(fn(f)  and  alf  =  -(s^nf  -s(fn[a>) . 


(2.4) 


2  '  2  ' 

As  a  measure  of  plastic  strain,  the  effective  plastic  shear  slip  is  calculated  from  the  plastic 

deformation  rate  tensor  as 
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N—p — f- 
D>D« 


Yeff 


(2.5) 


The  stress  is  updated  using  the  Jaumann  stress  rate  corotational  with  the  lattice,  of 


as 

af*=  C  De ,  (2.6) 

ij  ijkl  kl 

where  Cjkl  is  the  fourth-order  isotropic  elastic  modulus  tensor  defined  by 

cm  =  M(Slk8ji  +  SjkSu)  +  AS..Skl  (2.7) 

The  Jaumann  stress  rate  is  related  to  the  material  stress  rate,  cr,7 ,  in  the  reference  coordinate 
system  as 

<S>,e  e  e 

&v=(7v  +Wikakj  +  Wjkak,  (2.8) 

Power  law  hardening  is  assumed,  relating  the  slip  rates  on  the  various  slip  systems  to 
the  resolved  shear  stress  as 


/ a)  =  0 


no  sum  on  a. 


(2.9) 


&  re/)&  ref ) 


where  0  is  the  reference  shear  strain-rate  which  corresponds  to  a  reference  shear  stress, 


r";/ ,  and  m  is  the  strain-rate  sensitivity  parameter, 


_  <71n  r(g> 

^  • (a) 


din 


(2.10) 


The  reference  shear  stress  includes  forest  hardening  associated  with  the  immobile 
dislocation  density,  which  accounts  for  the  interaction  between  slip  systems  through  the 
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coefficient,  a • ,  and  thennal  softening  through  the  negative  thermal  softening  exponent,  £ , 
and  reference  temperature  Tr  (293  K)  [34], 


j=\,nss 


T* 


(2.11) 


For  high  strain-rate  investigations  under  the  assumption  of  adiabatic  heating,  the  temperature 
is  updated  using 

T  =  <j  Dp ,  (2.12) 

ij  ij 

PCP 

where  %  is  the  fraction  of  plastic  work  converted  to  heat,  p  is  the  mass  density,  c p  is  the 
specific  heat  of  the  material,  and  cri;  is  the  deviatoric  stress. 


2.2  Dislocation  Density  Evolution 

It  will  be  assumed  that  the  total  dislocation  density  can  be  decomposed  additively  into 
mobile  and  immobile  components. 


pH  =  A-a)  +  p(a)  (2.13) 

During  an  increment  of  slip,  mobile  dislocations  may  be  generated,  immobile  dislocations 
may  be  annihilated,  or  junctions  may  be  formed  or  destroyed  coupling  the  mobile  and 
immobile  dislocation  densities,  leading  to  the  coupled  differential  equations  governing 
dislocation  density  evolution, 


_  («)  _  .  («) 

Pm  —  P generation  interaction¬ 


al 4) 


.  («)  .  (4 

oim  =  "  P interaction  +  annihilation  • 


(2.15) 
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The  dislocation  density  evolution  follows  the  formulation  of  Shanthraj  and  Zikry 


[35].  Dislocation  density  generation  is  related  to  the  distance,  yback ,  traversed  by  a 
dislocation  emitted  from  a  source  with  density,  psource,  both  related  to  the  spacing  of 
immobile  forest  obstacles,  i.e.,  Frank  Read  mechanism.  The  average  velocity  of  mobile 
dislocations,  v^a\  is  used  to  detennine  the  generation  rate  as 


(a) 

P  generation  —  psource 


v(a) 


(2.16) 


y  back 


The  Orowan  equation,  j^a)  =  allows  the  generation  rate  to  be  recast  as 

'  («)* 

/^generation  — 11 - ) P‘m  ,yia) .  (2.17) 

b(a)  (pte)+ 

where  ^  is  a  geometric  parameter  and  bi")  is  the  magnitude  of  the  Burger’s  vector  on  slip 
system  a. 

Dislocation  density  interaction  involves  the  immobilization  of  mobile  dislocation 
segments  due  to  junction  fonnation  that  occurs  when  dislocation  densities  on  different  slip 
systems  interact.  The  frequencies  of  interaction  between  mobile  dislocation  densities  on  slip 

system  a  and  mobile  and  immobile  dislocations  on  slip  system  j3  are  defined  as  p„  p„  'v^ 

and  Pm  Pimv^  •>  where  the  relative  velocity  between  the  slip  systems,  v  ("//) ,  is  defined  using 
the  Orowan  equation  as 

=  (2.18) 

Pm  b  P,n  b 
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The  length  of  the  formed  junction  is  assumed  to  be  proportional  to  the  spacing  of  immobile 
dislocations  as 


(2.19) 


and  only  a  fraction  of  these  junctions,  f0,  are  stable.  The  rates  of  immobilization  of 
dislocation  densities  on  slip  system  a  due  to  mobile  and  immobile  dislocations  on  slip 
system  P  are  therefore  foP^'  p,!!'  h  and  f0pln]  p\„  lcv(a) ,  respectively.  The  rate  of 


immobilization  of  mobile  dislocation  densities  on  slip  system  a  is  therefore 


(«)  =f  ' p( P)j  jtP+p(a)l 


p(P)l  j/“)  +P  l  jX/9  (P)  y(") 

m  epp  m  e7 tPf’  + ,(a) 

p  (  b  b  +  p  b' 


(2.20) 


The  addition  of  immobile  dislocation  densities  due  to  interactions  also  takes  into 
account  the  possibility  of  dislocation  interaction  forming  immobile  junctions.  Frank’s  rule  is 
used  to  detennine  energetically  favorable  interactions  for  immobile  junction  formation.  A 
dislocation  density  interaction  tensor  for  junction  fonnation  on  slip  system  a  due  to 


interaction  of  slip  systems  J3  and  y  is  defined  as 


nSr=8f 


%lif  pb[a)1  <  pb{p)1  +  pb(y)1  and  b(a)  =b(p)  +b(r) 


0  otherwise 


(2.21) 


The  rates  of  junction  formation  on  slip  system  a  are  then  defined  as 

Pr  ( P )  M  (Pr)  Pr  (  (P)  (r)  (P)  (r)  (P)  W'j 

na  foPm  Pm  lcv  and  na  fpJc  Pm  pim  v  v  ' 


(2.22) 
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for  mobile/mobile  and  mobile/immobile  interactions.  The  total  addition  of  immobile 
dislocation  densities  due  to  interactions  then  becomes 
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Pi 


,(«) 


interaction  + 


+  Am  K 


( 


P,r 


b(r) 


2.23) 

py  b(r)  za>* 


* 


The  interaction  tensor  for  RDX  crystals  is  determined  using  the  6  energetically 
favorable  interactions:  self- interaction,  co-linear  interaction,  co-planar  interaction,  Lomer 
locks,  glissile  junctions,  Hirth  locks  [36].  For  RDX  crystals,  only  self-interaction  occurs. 
Table  2.1  summarizes  self-interaction  where  ay  is  assumed  to  be  0.6  only  for  self  interaction. 
Dislocation  density  annihilation  due  to  recovery  is  modeled  using  an  Arrhenius  relationship 
[37]  as 


Ml 

'  kT  * 


■Xa)  =  Me kr*  (2  24) 

H  annihilation  v  K 

where  the  frequency  of  which  immobile  dislocations  are  intersected  by  mobile  dislocations 
on  other  slip  systems  is  related  to  the  attempt  frequency,  ,  as 


p  c  b{p) 


(2.25) 


The  activation  enthalpy,  H,  is  related  to  the  immobile  dislocation  density  and  saturation 
density,  ps,  as 


h=h;  i- 


p\a)  ** 


(2.26) 


A  ) 


The  annihilation  rate  of  dislocation  densities  on  slip  system  a  becomes 


* 


>-#  V  K 

i 

kT 


.  («) 

P annihilation  j  0  /  . 

(  P 


pM  \  J 

im  c  ,  (t0  9 


(2.27) 
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The  generation,  interaction,  and  annihilation  rates  are  then  substituted  into  equations 
2.14  and  2. 1 5  to  obtain  a  coupled  nonlinear  set  of  evolutionary  equations  for  the  dislocation 


densities.  The  evolutionary  equations  can  be  expressed  as 

dp  0//°s  0//°  P  ( 

Hm  _  :/x  Ssour  Him 


_ y(x  <S  sour  H  im  v § immob—  P" 


a  a 


dt 

&b2 

cc  mnter- 

A 

&  m) 

dPim 

bo 

ii 

p’t  _|_  Simmob+ 

a 

a 

■i 


* 


a  a  a 

P  ~S  P 

'  im  ^  resow' 


dt 


mnter+  m 


n+ 


(2.28) 


(2.29) 


to  delineate  the  dislocation  activities  such  as  generation,  interaction,  immobilization,  and 
annihilation[33].  The  coefficients,  not  known  a  priori,  are  summarized  in  Table  2.2.  The 
dislocation  activity  is  coupled  to  the  stress  response  through  the  Taylor  relationship  (eq. 
2.11). 


2.3  Finite  Viscoelasticity  Formulation 

The  finite  viscoelasticity  formulation  for  the  polymer  binder  used  is  based  on  the 
approach  of  Kaliske  and  Rothert  [38],  Finite  viscoelasticity  is  employed  to  account  for  the 
general  viscous  nonlinear  behavior  of  estane  at  operating  temperatures  that  are  generally 
above  the  glass  transition  temperature.  The  glass  transition  temperature  of  estane  is 
approximately  -40°C  [39],  which  indicates  that  the  binder  would  be  in  a  viscoelastic  response 
state  for  the  prescribed  dynamic  loading  conditions  that  have  been  considered  in  this 
investigation.  The  finite  viscoelastic  constitutive  model  is  based  on  a  generalized  Maxwell 
model  (Figure  2.1).  It  is  assumed  that  the  stresses  can  be  decomposed  into  deviatoric  and 
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hydrostatic  components  as  the  deviatoric  component  exhibits  rate-sensitivity  and  the 
hydrostatic  component  exhibits  rate-insensitivity  as 
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(2.30) 


:  tt 

where  the  superscript,  n,  is  the  current  timestep. 

The  hydrostatic  stresses  are  updated  using  finite  elasticity  by  using  the  defonnation- 

rate  tensor  and  the  Jaumann  rate  of  Cauchy  stress  as 

(2.3T)  <C!"#  ■  Dr#  +  o'  ■ W-W-o-  , 

cr!"#  !"#  !"# 
where  Dl  #  is  the  hydrostatic  component  of  the  rate  of  deformation  tensor,  d'  "u  is  the 


hydrostatic  stresses  at  the  current  timestep,  W  is  the  total  spin  tensor,  and  (C!  u  is  the 
isotropic  elastic  material  matrix,  and  is  given  as 

rl  —  v  v  0 

(C!"#  =  !!  v  1-v  0 


(2.32) 


0  0  1  -  iV 


where  K  is  the  bulk  modulus,  and  v  is  Poisson’s  ratio. 

The  deviatoric  component  of  stress  is  updated  by  summing  the  contributions  from 

each  of  the  elements  in  the  Maxwell  model  as 

,,(2.33) 

where  cr,  r#  is  the  deviatoric  component  of  the  lone  spring  in  the  Maxwell  model,  H,  is 

the  viscous  stress  in  each  of  the  spring-dashpot  elements,  and  j  is  the  corresponding  number 
associated  with  each  spring-dashpot  in  the  Maxwell  model.  The  viscous  stresses  are  updated 

as 


=  <J—  + 


!"#  \\”  #  III 


V  V  V  V 

(  '  ), 

C\  i " i  cr |  j"# 


(2.34) 
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where  is  the  viscous  stress  at  the  current  timestep,  At  is  the  timestep  size,  T\  is  the 
H' 

relaxation  time,  G\  is  the  shear  modulus  of  the  spring-dashpot  element,  G\  is  the  shear 

!  !"# 

modulus  of  the  lone  spring,  and  o'  is  the  deviatoric  component  of  the  lone  spring  at  the 
current  timestep.  The  deviatoric  stresses  in  the  lone  spring  are  updated  using  a  finite 

elasticity  formulation  as  used  for  the  hydrostatic  component  (equations  2.30-2.32).  Time- 

temperature  superposition  is  used  to  modify  the  relaxation  times  through  the  use  of  a  shift 


factor,  di,  as 


log  (a 


(  )’ 

!"#!  !  ! !  i 


(2.35) 


where  T  is  a  reference  temperature  and  T  is  the  current  temperature. 


2.4  Thermo-Mechanical  Coupling 


The  heat  generation  due  to  mechanical  energy  is  given  as 

Dp, 

a  =y  (j .. 

*  mechanical  A  lJ  lJ 


(2.36) 


where  %  is  the  fraction  of  plastic  work  converted  to  heat,  and  o'  is  the  deviatoric  stress,  and 
D?j  is  the  plastic  deformation  rate  tensor.  The  heat  generated  due  to  thennal  decomposition  in 

the  energetic  regions  can  be  given  in  terms  of  an  Arrhenius,  exothennic  process  [40]  as 


-Ea 

q,  =PAHZeRT,  (2.37) 

1  decomposition  r  v  7 

where  p  is  the  mass  density,  AH  is  the  heat  of  decomposition,  Z  is  a  pre-exponential  factor, 

E,  is  the  energy  of  thermal  decomposition  activation,  R  is  the  universal  gas  constant,  and  T 
is  the  current  temperature.  The  heat  generated  due  to  viscous  dissipation  in  the  viscoelastic 

regions  can  be  given  by  summing  the  viscous  losses  over  the  entire  Maxwell  model  [41]  as 
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^  viscous  -  -  -  !  !  i  !  i 

For  high  strain-rate  simulations,  equations  (2.36  -  2.38)  can  then  be  used  to  detennine 


the  adiabatic  temperature  rate  as  follows, 


T  = 


l  !  i 


+ 


^  decomposition 


(2.39) 


where  c,  is  the  specific  heat.  Equations  (2.36  -  2.38)  can  also  be  used  as  heat  sources  to  the 
heat  conduction  equation  (2.40).  The  mechanical  energy  is  coupled  to  the  thennal  energy  and 


the  three  internal  heat  generation  sources  as  follows, 


PCp 


T+(in 


mechanical  ^ decomposition 


+ 


(2.40) 


2.5  Tables  and  Figures 


Table  2. 1 :  Summary  of  interactions  and  coefficient  values  in  an  RDX  crystal 


Interaction  type 

Dissipation  (oc yaij) 

ciij 

Self 

2kGb 7 

0.6 
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Table  2.2:  Summary  of  g-coefficients  for  dislocation  density  evolution  equations 


g  Coefficients 

Expression 
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CHAPTER  3:  Numerical  Methods 


The  stress  state  of  the  crystalline  materials  can  be  updated  by  means  of  the 
constitutive  relationship  prescribed  in  Chapter  2  once  the  total  velocity  gradient,  Ltj,  and 

plastic  velocity  gradient,  L1’.,  are  known.  The  specialized  finite  element  scheme 
implemented  for  the  detennination  of  the  velocity  gradient  tensors  will  be  outlined  below. 

3.1  Determination  of  the  Total  Velocity  Gradient 

The  total  velocity  gradient  can  be  obtained  by  means  of  finite  element  analysis.  For 

quasi-static  analysis,  an  incremental,  iterative  approach  taking  advantage  of  the  quasi- 
Newton  BFGS  scheme  is  used  to  detennine  the  nodal  displacements.  For  dynamic  analyses, 
an  implicit  Newmark-  j3,  iterative  approach  using  BFGS  to  solve  the  finite  element  equations 
linearized  about  t„  at  each  timestep  is  used.  Trapezoidal  values  of  /?  =  1/4  and  ;/=  1  [2  were 
chosen  for  unconditional  stability  of  the  finite  element  method.  The  deformation  calculated 
by  the  finite  element  method  is  used  to  calculate  the  total  velocity  gradient  in  accordance 
with  equation  2.1. 

To  avoid  numerical  locking  due  to  incompressible  pressure  constraints,  1 -point 
integration  of  the  Q4  quadrilateral  element  is  used,  which  has  the  added  benefit  of  reduced 
computational  time.  However,  reduced  integration  can  lead  to  the  zero-energy  numerical 
instability  of  hourglassing.  Stiffness-based  hourglass  control  is  implemented  to  control  the 
hourglass  instability. 
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Following  the  method  of  Flanagan  and  Belytschko  [42],  an  hourglass  base  vector 
hT=[l  -1  1  -1]  is  used  to  describe  the  hourglass  mode  for  a  one  point  integrated  Q4  element. 
The  hourglass  shape  vector,  which  is  orthogonal  to  the  linear  velocity  field  when  defined  in 
the  corotational  frame,  is 

h  =  ^[h-(hrx,)b],  i  =  x,y,  (3.1) 

where  b,  are  the  shape  function  derivatives,  and  x,  are  the  nodal  coordinates.  The 
generalized  hourglass  velocity  strains  are  given  in  terms  of  the  nodal  velocities,  v, ,  as 
(expressed  in  the  corotational  frame) 

and  qy  =  y'  v, .  (3.2) 

For  stiffness  based  hourglass  control,  generalized  hourglass  stress  rates  are  defined  as 

Q,  =  CQqx  and  Qy  =  C%,  (3.3) 

The  stress  rates  are  integrated  and  used  to  augment  the  force  vector  through  a  stabilization 
force, 

,  .  #Qxy& 

rtab=A3x  '  (3.4) 

%Qj( 

where  A  is  the  element  area.  A  corresponding  stabilization  augments  the  stiffness  matrix  as 

Ke=K  e,lpl  +  CQA!/  (3.5) 

3  yyT  yyT' 

The  coefficient  for  the  stabilization  force  and  stiffness  is  defined  as 

CQ  =  \ac2pAhTb  ,  (3.6) 

s  i  i 


22 


where  the  input  parameter  is  as,  and  c  is  the  shear  wavespeed  of  the  material.  An  effective 
modulus  is  used  to  account  for  the  lower  wavespeed  of  plastic  elements  as  follows: 


JlD 

<S>c t„/3.  -  LSt  ,7. 
D„dt  ’  H  2  DjDj  ’ 


A  =  K—  _ju 
3 


c 


(3.7) 


(3.8) 


3.2  Determination  of  the  Plastic  Velocity  Gradient 

The  objective  stress  rate  is  coupled  with  the  time  derivative  of  the  resolved  shear 

stresses  to  detennine  the  resolved  shear  stresses  on  each  slip  plane. 


(3.9) 

T'U)  =  C  p{a)f\e I  (3.10) 

ijkl 1  ij  ^ kl 


Eq.  3.10  can  be  expanded  to 


ir(a)  = 


2<P0Df-Z/’1 


ij  ref 


g=l,nss 


il 

1*®3 


($)\ 

Tref  + 


(3-11) 


where  the  reference  stress  is  a  function  of  the  immobile  dislocation  density.  The  slip  rates 
can  be  determined  using  the  resolved  shear  stress  (equation  2.9).  The  plastic  deformation 
rate  tensor  and  spin  rate  tensor  can  then  be  used  along  with  the  Schmid  tensor  and  slip  rates 
to  detennine  the  plastic  velocity  gradient  (equation  2.3). 
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The  elastic  spin  rate  tensor  can  be  detennined  as  a  function  of  the  total  spin  rate 
tensor  and  the  plastic  spin  rate  tensor  using  equation  2.2.  It  was  assumed  in  the  derivation  of 
equation  3.11  that  the  time  rate  of  change  of  the  slip  normals  and  directions  is  related  to  the 
lattice  spin  as 

«-a)  =  Wy'rij  and  sff>=W*Sj.  (3.12) 


3.3  Solving  Systems  of  ODE’s 

The  non-linear  ODE’s  governing  the  resolved  shear  stresses  and  dislocation  density 
evolution  necessitate  a  solution  method  that  is  both  accurate  and  stable  due  to  their  non-linear 
nature  and  the  possibility  of  numerical  stiffness  [43].  The  non-linear  nature  of  the  problem 
requires  a  solution  method  with  high  order  accuracy,  while  the  possibility  of  numerical 
stiffness  requires  a  solution  method  that  is  stable  and  will  not  propagate  error  due  to  stiff 
behavior. 

An  adaptive  timestep  fifth-order  accurate  step  halving  Runge-Kutta  method  is 
employed.  Two  approximate  solutions  are  taken  using  the  fourth-order  Runge-Kutta  method, 
one  at  step  length  h  and  a  second  with  two  steps  of  length  hf 2  . 

(6) 

r(t+h)=  zx+(ji)  <j)+0  h  ■+  (3.13) 


z(t+2(h  /2))  =  t2  +2(h  12)  (j)+ O 
The  two  solutions  are  combined  to  yield 


O 

h  +  *  ' 


(3.14) 


<S>i 


n 


z(t+h)  =  z2  +  +0  h 


(3.15) 
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where  <S>  ,  =  r,  —  f2  is  the  local  truncation  error,  which  is  used  to  measure  the  accuracy  of 
the  solution.  If  the  accuracy  is  not  less  than  a  specified  tolerance,  CED  0,  the  timestep  size 
is  reduced  by 


0.20 


/?,,  =  Fh 


old 


C>3 


(3.16) 


where  hnew  is  the  new  timestep  size,  hold  is  the  initial  timestep  size,  and  F  is  used  to  keep 
the  new  timestep  small  enough  to  be  accepted. 

When  the  timestep  is  becomes  excessively  small,  the  timestep  restriction  can  be  due 
to  stability  issues  caused  by  numerical  stiffness  [43].  The  integration  method  is  switched  to 


the  first-order  accurate,  unconditionally  stable,  backward  Euler  method, 

hf\^n+\Jn+ 1/5 
(a)  (a)  v  (a)  / 


(3.17) 


which  is  solved  using  quasi-Newton  iteration. 


3.4  Finite  Element  Representation  of  Thermal  Conduction 

The  discretized  finite  element  heat  conduction  equation  is  given  as 


[  ][']  [  ][  ]  [  ] 

C  T  +  K  T  =  R,  (3.18) 


[  ]  [t]  is  the 

where  C  is  the  matrix  of  rate  of  change  of  temperature  proportional  coefficients, 

vector  of  change  in  nodal  temperatures,  ^  is  a  matrix  of  temperature  proportional 


[  ]  [  ] 

coefficients,  T  is  the  vector  of  nodal  temperatures,  and  R  is  the  vector  of  nodal  input  heat 
sources  for  mechanical  and  decomposition  energies. 


The  elemental  matrix  of  rate  of  change  of  temperature  proportional  coefficients  is 


given  as 
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(3.19) 


1  =  /  U-  1W 


p  Ci  dV, 


[  ] 

where  N  is  the  vector  of  shape  functions,  and  the  integral  is  over  the  element  volume,  V. 
The  elemental  matrix  of  temperature  proportional  coefficients  is  given  as 


isotrop  jc  yi, 


aterial 


[  ] 


K 

[  ] 


k  0 

0  k 


The  vector  of  nodal  input  heat  sources  is  given  as 


V 


iCgtal  is  over  ,th  c-clcmen  tafi 
[  re  given  as  the  following  n 


m, 


owing  ior  an 


e,  V. 


(3.21) 


w 

VO 

V. 


llumc?;  's  ^1C  elemental  heat  get|ieratipn  rjt't(|/^in|l^hc integral  is  over  the  element!? 


The  temperatures  at  the  next  time  step  are  determined  using  the  following  equation, 

1  1 

(-[  ]  [  ]){  1  (-[  ]  (  )[  ])t  ! 
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where  At  is  the  time  step  size,  /?  was  taken  as  0.5  for  Crank-Nicolson  (trapezoidal) 
integration. 
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CHAPTER  4:The  Effects  of  Microstructural  Defects  on  Hot  Spot  Formation  in  RDX- 
PCTFE  Energetic  Aggregates 

4.1  Introduction 

When  the  thermo-mechanical  sensitivity  of  a  general  chemically  unstable  energetic 
material  is  large,  detonations  can  occur.  Hence,  understanding  and  controlling  the  detonation 
sensitivity  of  energetic  materials  is  critical.  It  is,  therefore,  essential  to  predict  and  identify 
inherent  dominant  microstructural  sources  of  sensitivity  due  to  thermo-mechanical  stimuli 
and  hot  spot  mechanisms. 

The  objective  of  this  investigation  is  to  predict  what  the  interrelated  microstructural 
effects  of  crystal-binder  interactions,  void  distributions,  dislocation-density  evolution,  and 
GBs  have  on  hot  spot  formation  in  aggregates  subjected  to  high  pressure  loading  conditions. 
Hence,  a  recently  developed  dislocation-density  based  crystalline  plasticity  approach  [35] 
and  specialized  dynamic  finite-element  techniques  will  be  used  to  predict  and  understand 
how  dislocation-density,  thermal  evolution,  inelastic  slip  behavior,  and  high  dynamic 
pressures  affect  hot  spot  fonnation  in  RDX-PCTFE  aggregates. 

4.2  Results 

Computational  analyses  were  undertaken  to  investigate  and  predict  the  effects  of 
different  dynamic  pressure  loading  conditions  and  GB  misorientations  on  hot-spot  fonnation 
for  a  RDX  energetic  crystalline  aggregate  with  a  polymer  binder.  The  finite  element 
aggregate  consisted  of  25  grains  of  RDX  crystals  with  an  average  grain  size  of  200  pm  for  an 
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aggregate  size  of  1mm  x  1mm  (Figure  4.1).  The  RDX  crystalline  morphology  was  generated 
utilizing  a  Voronoi  tessellation  scheme.  Based  on  Gallagher,  Halfpenny,  et  al.  [44],  the  RDX 
grains  were  assumed  to  have  the  following  three  slip  systems  -  (010)  [001],  (021)  [100], 

(02l)  [100]  -  with  material  properties  given  by  Armstrong  and  Elban  [1].  A  polymer  binder 
of  Polychlorotrifluoroethylene  (PCTFE),  with  a  volume  fraction  of  10%,  was  randomly 
distributed  between  the  RDX  grains.  The  polymer  properties,  which  are  based  on 
measurements  by  Brown,  Rae,  et  al.  [45],  are  given  in  Table  II.  The  polymer  binder  has  been 
modeled  as  a  purely  elastic  material  with  a  high  strength  to  ensure  plasticity  does  not  occur 
because  the  glass  transition  temperature  is  approximately  50°C  [45],  which  indicates  that  the 
binder  would  be  in  the  glassy  state  for  the  dynamic  loading  conditions  that  have  been 
considered  in  this  investigation. 

Porosity  in  RDX  crystals  can  be  controlled  by  changing  the  crystal  processing 
parameters,  which  yields  many  different  combinations  of  void  sizes,  distributions,  and 
overall  porosity  volume  fractions  that  are  possible  [7,  46].  In  this  investigation,  circular 
voids,  with  an  average  diameter  of  15  pm  and  an  overall  porosity  volume  fraction  of  3% 
were  used  (Figure  4. 1).  Symmetry  boundary  conditions  were  applied  to  the  model  for  two- 
dimensional  plane-strain  loading  conditions.  A  constant  dynamic  pressure  load  was  applied 
on  the  top  surface,  as  shown  in  Figure  4.1,  for  pressure  loads  that  varied  from  750  MPa  to  1 
GPa  for  random  low  (misorientations  less  than  15°)  and  random  high  angle  GB 
misorientations  (misorientations  greater  than  30°). 
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4.2.1  750  MPA  Pressure  Loading  and  Random  Low  Angle  GB  Distributions 

For  a  compressive  pressure  of  750  MPa,  at  a  time  of  630  nanoseconds,  with  random 

low  angle  GB  misorientations,  the  lattice  rotation,  which  is  a  measure  of  geometrical 
softening,  varied  from  +60°  to  -60°  (Figure  4.2(a))  where  the  maximum  values  occurred  at 
the  periphery  of  the  voids.  The  normal  stress  was  generally  compressive  in  the  bulk  of  the 
material.  There  were,  however,  localized  regions  at  the  top  and  bottom  surfaces  of  the  void 
peripheries  where  tensile  nonnal  stresses  occurred  (Figure  4.2(b)).  The  stresses  and  pressure 
were  normalized  by  the  static  yield  strength  of  RDX.  The  pressure  was  nominally 
compressive  in  the  bulk  of  the  material  with  regions  around  the  entire  void  periphery  that 
were  tensile  (Figure  2(c)).  The  tensile  normal  stresses  and  pressures  most  likely  occured  due 
to  wave  reflections  and  interactions  around  the  void  peripheries  [47].  The  accumulated 
plastic  shear  strain  (Figure  2(d)),  which  is  the  total  inelastic  slip  on  all  active  systems,  had 
also  localized  around  the  peripheries  of  the  voids  with  a  maximum  value  of  160%. 

The  immobile  dislocation  densities  for  the  three  independent  slip  systems  in  RDX  are 
shown  in  Figure  4.3(a-c).  The  dislocation  density  for  the  (010)  [001]  slip  system,  slip  system 
one  in  Figure  4.3(a),  had  a  maximum  nonnalized  value  of  350.  The  dislocation  densities 
were  normalized  by  the  initial  dislocation  density.  The  dislocation  densities  for  slip  systems 

(021)  [100],  slip  system  two,  and  (021)  [100],  slip-system  three  in  Figure  4.3(b-c),  have 
maximum  normalized  values  of  30.  This  indicates  that  for  random  low  GB  misorientations, 

slip  system  one  was  clearly  dominant,  and  plastic  deformation  occurred  by  dislocation  slip 

primarily  due  to  this  slip  system. 
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Figure  4.4(a)  shows  the  spatial  temperature  distribution,  which  includes  the  effects  of 
thennal  conduction,  plastic  work,  and  thennal  decomposition.  The  temperature  accumulated 
in  the  intervoid  regions  where  multiple  void  deformations  interacted  and  where  the  shear 
strains  have  localized,  and  the  maximum  temperature  in  this  region  was  460K.  Therefore,  the 
maximum  temperature  rise  in  the  aggregate  was  167K  (the  reference  temperature  was  293K) 
in  the  intervoid  regions.  Figure  4.4(b)  shows  the  temperature  rise  profile  due  to 
decomposition,  which  indicates  that  thennal  decomposition  had  occurred  for  this  pressure 
loading  and  GB  distribution.  The  maximum  decomposition  temperature  rise  was  3K,  which 
occurred  at  the  periphery  of  one  specific  void,  and  the  maximum  adiabatic  plastic  work 
temperature  (Figure  4.4(c))  was  700K.  This  substantiates  that  even  at  small  time  scales  and 
high  loading  rates,  the  effects  of  thermal  conduction  were  not  negligible,  because  the 
maximum  adiabatic  plastic  work  temperature  (Figure  4.4(c))  was  not  comparable  to  the 
maximum  temperature  of  the  aggregate  (Figure  4.4(a))  that  included  heat  conduction.  If  these 
two  temperatures  were  of  comparable  magnitude,  then  the  assumption  of  thermal  conduction 
being  negligible  would  have  been  physically  acceptable. 

The  temperature  and  normal  stress  evolution,  as  a  function  of  time,  for  the  element 
that  had  the  highest  temperature  rise  is  shown  in  Figure  4.5.  This  element  was  between  two 
voids,  where  the  shear  strain  had  localized.  As  the  plot  indicates,  the  temperature  was 
constant  until  plastic  deformation  occurred,  and  then  it  increased  sharply  to  460K  over  a 
time-span  of  approximately  20  nanoseconds.  The  nonnal  stress  in  this  region  attained  a 
maximum  normalized  compressive  value  of  -7,  the  temperature  increased,  and  the  normal 
stress  increased  rapidly  to  a  compressive  value  of  -2.  According  to  Bruckman  and  Guillet[48] 
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and  Field,  Borne,  et  al.  [5]  a  hot  spot  can  form  or  initiate  when  the  heat  generated  due  to 
decomposition  exceeds  the  heat  dissipated  to  the  surroundings,  thus  causing  a  run-away 
temperature  build-up.  Therefore,  based  on  these  results,  we  can  predict  hot  spot  initiation 
when  the  local  temperature  becomes  unbounded  due  to  decomposition.  For  these  loading 
conditions  and  GB  misorientations,  the  temperature  (Figure  4.5)  was  bounded,  and  therefore 
it  is  predicted  that  hot  spot  initiation  would  not  occur. 

4.2.2  1  GPA  Pressure  Loading  and  Random  Low  Angle  GB  Distributions 

For  a  compressive  loading  of  1  GPa,  at  a  time  of  430  nanoseconds,  with  random  low 

angle  misorientations,  the  lattice  rotations  (Figure  4.6(a))  varied  between  +60°  and  -60°, 
where  the  maximum  values  occurred  around  the  void  peripheries.  The  normal  stress  (Figure 
4.6(b))  had  a  maximum  nonnalized  tensile  stress  of  3  and  a  normalized  maximum 
compressive  stress  of  8.  The  aggregate  was  nominally  compressive,  but  tensile  stress  regions 
occurred  at  the  top  and  bottom  surfaces  of  the  void  peripheries,  which,  as  in  the  750  MPa 
case,  was  likely  due  to  wave  reflections  and  wave  interactions  around  void  peripheries. 
Figure  4.6(c)  shows  the  pressure  distribution  throughout  the  aggregate.  The  maximum 
nonnalized  tensile  pressure  was  1.5,  and  the  maximum  nonnalized  compressive  pressure  was 
4.5.  The  maximum  accumulated  plastic  shear  strain  was  100%,  and  it  occurred  adjacent  and 
around  the  void  peripheries  (Figure  4.6(d)),  and  it  is  less  than  the  maximum  shear  strain  for 
the  750  MPa  loading  by  approximately  60%.  This  decrease  in  plastic  slip  was  likely  due  to 
the  increased  strain  hardening  associated  with  the  higher  dynamic  loading. 

The  immobile  dislocation  density  for  slip  system  one  (Figure  4.7(a)),  (010)  [001], 
attained  a  maximum  nonnalized  value  of  160,  while  the  dislocation  densities  for  slip  systems 
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two  (Figure  4.7(b)),  (021)  [100],  and  three  (Figure  4.7(c)),  (021)  [100],  had  maximum 
nonnalized  values  of  16,  and  the  maximum  dislocation  densities  accumulated  around  the 

void  peripheries.  Slip  system  one  is  clearly  dominant,  which  as  with  the  lower  pressure  case 

of  750  MPa,  indicates  that  this  slip-system  clearly  dominated  the  overall  plastic  defonnation. 

The  temperature  (Figure  4.8(a))  had  a  maximum  value  of  520K,  the  maximum 
decomposition  temperature  rise  (Figure  4.8(b)),  before  becoming  unbounded,  was  13K,  and 
the  temperature  due  to  plastic  work  (Figure  4.8(c))  was  580K.  The  temperature  temporal 
evolution  (Figure  4.9)  in  the  element  where  the  temperature  was  the  largest  remained 
constant  until  plastic  deformation  occurred  at  which  point  the  temperature  increased  sharply 
to  520K.  Then  the  temperature  became  unbounded  due  to  thennal  decomposition,  and  this 
indicated  hot  spot  initiation.  The  temperature  became  unbounded  450  nanoseconds  after  the 
application  of  the  compressive  pressure  load.  The  normal  stress  was  initially  zero,  until  the 
stress  wave  reached  the  critical  element,  at  which  point  the  stress  increased  to  a  normalized 
compressive  value  of  1.5.  When  the  temperature  increased  due  to  plastic  defonnation,  the 
stress  transitioned  from  compressive  to  tensile.  The  hot  spot  formation,  at  this  higher 
pressure  loading  case,  is  in  contrast  with  the  750  MPa  loading  case  (Figure  4.5),  where  the 
temperature  was  bounded,  and  therefore  hot  spot  initiation  would  not  have  occurred.  This 
would  further  indicate  that  as  thennal  decomposition  becomes  more  dominant  with  increases 
in  dynamic  pressure  loading,  hot  spot  initiation  is  more  likely  to  occur. 
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4.2.3  750  MPA  Pressure  Loading  and  Random  High  Angle  GB  Distributions 

For  a  compressive  loading  of  750  MPa,  at  a  time  of  600  nanoseconds,  with  random 

high  angle  GB  misorientations,  the  lattice  rotation  (Figure  4.10(a))  was  bounded  by  +60°  and 
-60°,  similar  to  the  low  GB  misorientations  with  a  pressure  load  of  750  MPa.  The  normal 
stresses  and  pressures  (Figure  4.10(b-c))  were  generally  compressive  in  the  bulk  of  the 
aggregate  with  localized  regions  of  tensile  stresses  that  occurred  around  the  void  peripheries. 
The  accumulated  plastic  shear  strain  (Figure  4. 10(d))  had  a  maximum  of  130%  within 
localized  deformation  regions. 

The  dislocation  density  for  slip  system  one  (Figure  4. 1 1(a)),  (010)  [001]  had  a 
maximum  normalized  value  of  240.  The  dislocation  density  for  slip  system  two  (Figure 
4.1 1(b)),  (021)  [100],  had  a  maximum  normalized  value  of  100,  and  slip  system  three  (Figure 

4.1 1(c)),  (021)  [100],  had  a  maximum  normalized  value  of  26.  Slip  system  one  is  clearly  the 
dominant  slip  system,  but  slip  system  two  had  more  activity  in  comparison  with  the  random 

low  angle  GB  case  at  a  loading  of  750  MPa,  which  is  an  indication  that  slip  behavior  is  no 

longer  dominated  solely  by  one  slip  system,  as  was  the  case  for  the  random  low  angle  GB 

cases.  The  activation  of  this  slip  system  can  be  an  indication  of  larger  stresses  and  stress 

gradients  that  cross  between  the  binder  and  the  RDX  crystals,  and  these  higher  stresses 

would  result  in  the  activation  of  slip  systems. 

The  maximum  temperature  of  550K  (Figure  4. 12(a))  occurred  at  the  intervoid  region 

as  indicated  in  the  figure.  The  maximum  decomposition  temperature  rise  (Figure  4.12(b)) 

was  77K,  which  occurred  in  the  localized  defonnation  region.  The  maximum  temperature 

rise  due  to  plastic  work  (Figure  4. 12(c))  was  357K,  which  also  occurred  in  the  intervoid 
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localized  deformation  region.  Based  on  these  results,  it  is  clear  that  GB  misorientations 
affected  hot  spot  formation  since  hot  spot  formation  did  not  occur  with  lower  GB 
misorientations  for  the  same  loading  conditions.  This  is  likely  due  to  the  possibility  that  slip 
can  be  impeded,  and  that  high  stresses  would  then  occur. 

The  temporal  evolution  of  the  temperature  (Figure  4.13)  in  the  critical  element  where 
the  temperature  rise  was  the  largest  was  initially  constant  for  450  nanoseconds  at  which  time 
the  temperature  gradually  increased  until  it  became  unbounded  approximately  600 
nanoseconds  after  the  application  of  loading.  In  conjunction  with  the  thermal  buildup,  the 
nonnal  stress  was  initially  zero  until  200  nanoseconds  at  which  point  it  increased 
compressively  to  a  normalized  value  of  2.  As  the  temperature  increased,  the  normal  stress 
transitioned  from  compressive  to  tensile.  It  is  important  to  note  that  the  temperature  became 
unbounded  600  nanoseconds  after  the  application  of  the  applied  load,  and  this  would  then 
result  in  hot  spot  initiation  in  the  critical  element. 

Compared  with  experimental  studies  [2,  3,  7,  8,  11,  49],  our  predictions  further 
substantiate  the  relation  between  void  population,  dynamic  pressure  loading  conditions,  and 
the  sensitivity  to  hot  spot  formation.  Hot  spot  formation  would  occur  more  rapidly  with 
increases  in  incident  pressure.  Our  predictions  have  indicated  that  at  lower  pressure  loads,  the 
temperature  was  bounded,  and  therefore,  hot  spot  initiation  would  not  occur.  Conversely,  for 
higher  dynamic  pressure  loads,  unbounded  temperatures  would  result  in  hot  spot  initiation, 
and  this  is  also  consistent  with  the  above  investigations. 
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4.3  Summary 

A  dislocation-density  based  crystalline  plasticity  and  specialized  finite-element 
formulations  were  used  to  investigate  hot  spot  formation  due  to  dynamic  thermo-mechanical 
loading  conditions  in  RDX-PTCFE  energetic  aggregates.  At  a  dynamic  pressure  load  of  750 
MPa,  for  the  aggregate  with  low  angle  GB  misorientations,  the  bulk  of  the  aggregate  was 
generally  compressive  with  localized  regions  of  tensile  pressure  around  the  void  peripheries. 
Accumulated  plastic  shear  strain  deformation  occurred  in  localized  regions  around  the  void 
peripheries,  and  this  plastic  deformation  was  mainly  due  to  the  dominance  by  one  slip- 
system.  The  temperature  also  accumulated  in  regions  around  the  voids  where  the  shear  strain 
had  localized.  For  this  pressure  loading  and  GB  misorientation  distribution,  hot  spot 
formation  did  not  occur  since  there  was  no  extensive  thennal  decomposition  and  no 
unbounded  increases  in  temperature. 

For  the  1  GPa  loading  with  low  GB  misorientations,  the  shear  strain  deformation  had 
also  localized  around  the  void  peripheries,  and  the  plastic  deformation  was  also  due  to  one 
dominant  slip-system.  The  temperature  built  up  in  intervoid  regions,  and  became  unbounded 
due  to  increases  in  the  thermal  decomposition.  These  unbounded  temperatures  would  result 
in  hot  spot  initiation. 

For  the  750  MPa  compressive  loading  with  random  high  angle  GB  misorientations 
the  shear  strain  localized  in  banded  regions,  and  plastic  deformation  occurred  primarily  on 
two  of  the  three  slip-systems.  The  temperature  accumulated  in  the  banded  defonnation 
regions  where  the  shear  strain  had  localized.  The  temperature  became  unbounded,  due  to 
thennal  decomposition,  which  would  also  result  in  hot  spot  initiation.  Plastic  deformation 
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occurred  by  the  domination  of  dislocation-density  slip  on  two  slip  systems,  whereas  for  the 
low  misorientation  simulation,  plastic  deformation  occurred  dominantly  on  the  first  slip 
system.  The  high  misorientation  simulation  also  indicates  that  hot  spot  initiation  will  occur, 
whereas  for  low  misorientations,  hot  spot  initiation  would  not  occur.  The  differences  in 
sensitivity  of  the  high  and  low  angle  GB  misorientation  cases  indicates  that  an  energetic 
material  that  has  its  loading  axis  oriented  with  more  available  slip  systems  could  be  used  to 
lower  the  sensitivity  of  the  energetic  material  to  hot  spot  fonnation.  The  high  angle  GB  case 
also  indicated  that  high  angle  GBs  may  impede  slip  on  a  specific  slip-system,  which  would 
increase  local  stresses,  which  can  then  activate  other  slip-systems. 

Results  from  this  study  indicate  that  the  interrelated  effects  of  dynamic  pressure 
loading,  GB  distributions,  and  void  populations  affect  thennal  decomposition  and  hot  spot 
initiation.  Unbounded  temperatures  resulted  from  heat  released  during  the  decomposition 
reaction,  which  was  thermally  activated  by  heat  generation  from  plastic  defonnation.  Hot 
spot  fonnation  resulted  because  the  heat  generated  by  plastic  defonnation  and  decomposition 
was  greater  than  the  heat  dissipated  to  the  surrounding  material  through  thennal  conduction. 
This  indicates  that  for  certain  dynamic  pressure  loading  conditions,  the  heat  generated  by 
plastic  deformation  can  provide  the  thennal  activation  required  for  hot  spot  initiation  by 
thennal  decomposition. 
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4.4  Tables  and  Figures 


Figure  4.1:  The  RDX -polymer  aggregate  with  loading,  boundary  conditions,  and  void 

distribution 
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Figure  4.2.  (a)  Rotation,  (b)  Normal  Stress,  (c)  Pressure,  and  (d)  Accumulated  plastic  shear 
strain  for  a  750  MPa  compressive  loading  with  low  misorientations. 


38 


(C) 


Dislocation  Density 
Immobile  3 


0.0012 


Figure  4.3.  (a)  Immobile  dislocation  density  l,  (b)  Immobile  dislocation  density  2,  (c) 
Immobile  dislocation  density  3  for  a  750  MPa  compressive  loading  with  low  misorientations. 
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Figure  4.4.  (a)  Total  temperature,  (b)  Decomposition  temperature,  (c)  Plastic  work 
temperature  for  a  750  MPa  compressive  loading  with  low  misorientations. 
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Figure  4.5.  Temperature  and  normal  stress  as  a  function  of  time  for  a  750  MPa  compressive 

loading  with  low  GB  misorientations. 
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Figure  4.6.  (a)  Rotation,  (b)  Normal  Stress,  (c)  Pressure,  and  (d)  Accumulated  plastic  shear 
strain  for  a  1  GPa  compressive  loading  with  low  angle  GB  misorientations. 
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Figure  4.7.  (a)  Immobile  dislocation  density  1,  (b)  Immobile  dislocation  density  2,  (c) 
Immobile  dislocation  density  3  for  a  1  GPa  compressive  loading  with  low  GB 

misorientations. 
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Figure  4.8.  (a)  Total  Temperature,  (b)  Decomposition  Temperature,  (c)  Plastic  Work 
Temperature  for  a  1  GPa  Compressive  loading  with  low  misorientations. 
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Figure  4.9.  Temperature  and  normal  stress  as  a  function  of  time  for  a  1  GPa  compressive 
loading  with  low  angle  GB  misorientations.  Hot  spot  initiation  occurs  due  to  an  unbounded 
temperature  caused  by  increases  in  thermal  decomposition. 
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Figure  4.10.  (a)  Rotation,  (b)  Normal  Stress,  (c)  Pressure,  and  (d)  Accumulated  plastic  shear 
strain  for  a  750  MPa  compressive  loading  with  high  angle  GB  misorientations. 
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Figure  4.1 1.  (a)  Immobile  dislocation  Density  1,  (b)  Immobile  dislocation  density  2,  (c) 
Immobile  dislocation  density  3  for  a  750  MPa  compressive  loading  with  high  angle  GB 

misorientations. 
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Figure  4.12.  (a)  Total  Temperature,  (b)  Decomposition  Temperature,  (c)  Plastic  Work 
Temperature  for  a  750  MPa  Compressive  loading  with  high  misorientations. 
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Figure  4.13.  Temperature  and  normal  stress  as  a  function  of  time  for  a  750  MPa  compressive 

loading  with  high  misorientations. 
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CHAPTER  5:  Heterogeneous  Thermo-Mechanical  Behavior  and  Hot  Spot  Formation  in 
RDX-Estane  Energetic  Aggregates 

5.1  Introduction 

The  heterogeneous  crystalline-amorphous  structure  of  energetic  materials,  the  defects 
that  span  different  spatial  scales,  and  the  dynamic  thermo-mechanical  loading  conditions 
underscore  that  a  fundamental  understanding  of  behavior  and  hot  spot  formation  is  needed. 
Therefore,  the  objective  of  this  investigation  is  to  formulate  a  predictive  framework  that 
accounts  for  the  interrelated  microstructural  effects  of  crystalline  plasticity,  polymer  binder 
relaxation,  and  inherent  microstructural  defects  and  features  such  as  void  distribution, 
dislocation-density  evolution,  and  GBs  on  hot  spot  fonnation  mechanisms  in  aggregates 
subjected  to  high  pressure  mechanical  loads.  A  recently  developed  dislocation-density  based 
crystalline  plasticity [3 5]  approach,  which  has  been  coupled  to  a  finite  viscoelasticity 
formulation  and  specialized  dynamic  finite-element  techniques  will  be  used  to  predict  hot 
spot  fonnation  mechanisms  in  RDX-estane  aggregates. 

5.2  Results 

Computational  analyses  were  undertaken  to  investigate  and  predict  the  effects  of 
dynamic  loading  conditions  on  system  behavior  and  hot  spot  formation  for  a  RDX  energetic 
crystalline  aggregate  with  a  viscoelastic  polymer  binder.  The  finite  element  aggregate 
consisted  of  64  grains  of  RDX  crystals  with  an  average  grain  size  of  125  pm  for  an  overall 
aggregate  size  of  1mm  x  lnnn  (Figure  5.1).  The  RDX  crystalline  morphology  was  generated 
utilizing  a  Voronoi  tessellation  scheme.  Based  on  Gallagher  et  al.  [44],  the  RDX  grains  were 
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assumed  to  have  the  following  three  independent  slip  systems  -  (010)  [001],  (021)  [100], 
(02-1)  [100]  -  with  material  properties  given  by  Armstrong  and  Elban[l]  and  Annapragada  et 
al.  [50]  in  Table  5.1.  A  polymer  binder  of  estane,  with  volume  fraction  of  approximately 
17%,  was  randomly  distributed  between  the  RDX  grains.  The  material  properties  for  the 
estane  polymer  binder  were  obtained  from  measurements  by  Mas  et  al.  [39]  and  Barua  and 
Zhou  [23],  as  shown  in  Table  5.1. 

Porosity  in  RDX  crystals  can  be  controlled  by  changing  the  crystal  processing 
parameters,  which  can  result  in  different  combinations  of  void  morphology,  sizes,  and 
distributions  [7,  46].  In  this  investigation,  circular  voids,  with  an  average  diameter  of  10  pm 
and  an  overall  porosity  volume  fraction  of  2.5%  were  used  (Figure  5.1),  which  is  consistent 
with  energetic  material  microstructures  [7,  46,  49,  5 1,  52],  Symmetry  boundary  conditions 
were  applied  to  the  model  for  two-dimensional  plane-strain  loading  conditions.  A  constant 
dynamic  pressure  load  was  applied  on  the  top  surface,  as  shown  in  Figure  5. 1,  for  a  dynamic 
tensile  pressure  load  of  1  GPa  and  a  dynamic  compressive  pressure  load  of  1  GPa,  with 
random  GB  misorientations  of  less  than  15°. 

5.2.1  Validation  of  the  Viscoelastic  Model 

For  the  estane  polymer  binder,  the  model  was  validated  with  the  experimental  results 

of  Mas  et  al.  [39]  for  quasi-static  and  dynamic  behavior.  A  pure  estane  model  with 
dimensions  of  1mm  x  1mm  was  used,  and  symmetry  boundary  conditions  were  applied  for 
two-dimensional  plane-strain  loading  conditions.  A  constant  strain  rate  displacement 
condition  was  applied  on  the  top  surface,  as  shown  in  Figure  5.2,  for  a  quasi-static  strain  rate 
of  0.005  s'1  and  a  dynamic  strain  rate  of  2000  s'1. 
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For  a  quasi-static  strain  rate  of  0.005  s'1  and  an  initial  temperature  of  20°C,  predicted 
and  experimental  results  were  both  non-linear  (Figure  5.3a).  The  predicted  results,  at  a 
nominal  strain  of  15%,  were  within  approximately  10%  of  the  experimental  value,  and  the 
percent  difference  of  the  average  stresses  was  6.5%.  For  the  dynamic  strain  rate  of  2000  s'1 
and  an  initial  temperature  of  0°C,  the  percent  average  difference  between  the  experimental 
and  predicted  results  was  approximately  6.9%  (Figure  5.3b). 

The  differences  in  the  experimental  results  are  likely  due  to  the  experimental 
cylindrical  specimen  geometry,  while  the  analysis  was  for  a  plane  strain  simulation. 
Furthennore,  the  specimen  size  was  never  specified  in  [39],  therefore,  a  size  of  1  mm  x  1 
mm,  that  is  consistent  with  a  split  Hopkinson  specimen  was  used  (see,  for  example[53,  54]). 
The  viscoelastic  stress  evolutions  as  a  function  of  strain,  however,  were  consistent  with  both 
the  quasi-static  and  dynamic  experimental  results  and  the  average  differences  were  less  then 
10%.  Hence,  this  validation  is  another  indication  of  the  predictive  capabilities  of  the 
computational  approach. 

5.2.2  Dynamic  Tensile  Pressure  Loading  of  1  GPa 

Local  microstructural  behavior  was  investigated  by  using  a  64  grain  aggregate  that 

was  limn  x  1mm,  which  was  subjected  to  a  dynamic  tensile  pressure  load  of  1  GPa.  At  a 
time  of  .465  ps,  the  spatial  temperature  distribution,  which  included  the  effects  of  thermal 
conduction,  adiabatic  heating,  and  thennal  decomposition,  is  shown  in  Figure  5.4a.  The 
temperatures  accumulated  at  the  crystal  edges  where  the  shear  strains  had  localized,  as  well 
as  in  the  highly  sheared  regions  of  the  polymer  binder,  and  the  maximum  normalized 
temperature  in  this  region  was  1.87  (the  temperatures  were  nonnalized  by  the  initial 
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temperature  of  293  K).  The  maximum  temperature  rise  in  the  aggregate  was  257  K  over  a 
span  of  .465  ps.  The  normal  stress  was  generally  tensile  in  the  aggregate  (Figure  5.4b)  with  a 
maximum  normalized  stress  (the  stresses  were  normalized  by  the  yield  stress)  of  6.  However, 
there  were  regions  of  compressive  nonnal  stress  within  the  aggregate  that  occurred  at  the 
crystal  edges  and  in  the  adjoining  binder  with  a  maximum  normalized  compressive  normal 
stress  of  6.  The  compressive  normal  stress  at  the  crystal  edges  and  binder  occur  due  to  the 
viscosity  of  the  polymer,  which  facilitated  crystal  sliding  and  crystal  interaction  with  other 
crystals.  The  lateral  stress  was  generally  tensile  in  the  horizontal  ligaments  of  binder  (Figure 
5.4c)  with  a  maximum  value  of  4.  A  major  region  of  compressive  lateral  stress  occurred  at 
the  right  edge  of  the  model.  This  is  due  to  the  effects  of  the  free  boundary  and  associated 
release  waves.  High  regions  of  shear  stress  were  localized  at  the  crystal  edges  and  the 
associated  binder  regions  that  had  extensively  sheared  (Figure  5.4d).  The  maximum 
normalized  shear  stress  was  2.5  and  the  minimum  value  was  -2. 

Figure  5.5a  shows  the  spatial  distribution  of  dissipated  energy  due  to  plasticity  and 
viscous  effects  in  the  polymer  binder.  The  dissipated  energy  was  localized  at  the  crystal  edges 
where  the  temperatures  had  accumulated  with  a  maximum  value  of  1 .8E9  joules.  The 
hydrostatic  pressure  was  generally  tensile  in  the  aggregate  with  a  maximum  normalized  value 
of  3  occurring  in  the  binder  regions  that  were  fully  constrained  by  neighboring  RDX  grains 
(Figure  5.5b).  The  minimum  normalized  pressure  of  -6  occurred  in  the  binder  region  that  was 
under  compression  due  to  crystal  sliding.  These  regions  of  high  compressive  pressure  are  due 
to  the  dominance  of  hydrostatic  stress  in  the  binder  (see  Equations  10  and  11).  The  plastic 
deformation,  therefore,  at  the  neighboring  crystal  edges,  occurs  due  to  the  incompressibility 
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and  relative  ease  of  polymer  shearing.  This  underscores  how  crystal  binder  interfaces  play  an 
integral  role  in  the  aggregate  deformation.  Crystal  sliding  has  been  shown  to  be  important  in 
the  defonnation  of  crystalline  aggregates  as  predicted  by  Van  Swygenhoven  and  Derlet  [55] 
and  Kim  et  al.  [56]  and  this  results  in  stress  accumulations.  These  results,  furthermore,  are 
consistent  with  the  computational  modeling  results  of  Barua  et  al.  [22,  23],  which  indicate 
that  hot  spots  occur  due  to  intense  regions  of  shear  defonnation  in  the  binder  for  an  HMX 
and  estane  aggregate.  The  accumulated  plastic  shear  strain  (Figure  5.5c),  which  is  the  total 
inelastic  slip  on  all  slip  systems,  had  a  maximum  of  80%,  which  occurred  at  crystal  edges 
where  the  maximum  temperatures  accumulated.  This  results  in  heat  from  the  binder,  due  to 
viscous  dissipation,  diffusing  into  the  RDX  crystal  and  thus  increasing  the  temperature  and 
crystalline  plastic  defonnation  at  crystal-estane  interfaces. 

The  mobile  dislocation  densities  for  the  three  independent  slip  systems  in  RDX  are 
shown  in  Figure  5.6  (a-c).  The  dislocation  density  for  the  (010)  [001]  slip  system,  slip  system 
one  (Figure  5.6a),  had  a  maximum  nonnalized  value  of  8000  (the  dislocation  densities  were 
nonnalized  by  the  initial  dislocation  density),  which  occurred  at  the  crystal  edges  where  the 
temperatures  accumulated.  The  dislocation  density  for  slip  systems  (021)  [100],  slip  system 
two  (Figure  5.6b),  had  a  maximum  nonnalized  value  of  1000  and  the  dislocation  density  for 

(021)  [100],  slip-system  three  (Figure  5.6c),  had  a  maximum  normalized  value  of  22. 

The  immobile  dislocation  densities  are  shown  in  Figure  5.7(a-c).  The  dislocation 

density  for  the  slip  system  one  (Figure  5.7a)  had  a  maximum  normalized  value  of  800.  The 
dislocation  density  for  slip  system  two  (Figure  5.7b)  had  a  maximum  nonnalized  value  of 
100  and  the  dislocation  density  for  slip-system  three  (Figure  5.7c)  had  maximum  normalized 
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values  of  6.  This  indicates  that  for  low  angle  misorientations,  the  inelastic  deformation  is 
dominated  by  slip  system  one.  Experiments  by  Hooks  et  al.  [57]  and  Gallagher  et  al.  [44] 
indicate  that  RDX  deformation  is  highly  anisotropic  and  orientation  dependent.  This  further 
substantiates  that  for  random  low  GB  misorientations  that  one  slip  system  can  be  the 
dominant  slip  system.  Hence,  for  the  random  low  GB  misorientations  used  in  this  study,  the 
predictions  indicate  that  the  primary  pathway  for  plastic  defonnation  is  dominated  by  one 
slip  system,  and  based  on  the  experimental  observations  and  measurements,  this  is  due  to  the 
highly  orientation-dependent  dominance  of  specific  RDX  slip-systems. 

The  temperature,  normal  and  shear  stresses,  and  pressure  temporal  evolutions  for  the 
region  where  the  temperature  becomes  unbounded  is  shown  in  Figure  5.8.  The  spatial 
location  on  this  region  is  given  in  Figure  5.4a.  The  temperature  was  initially  constant  until 
plastic  deformation  initiated  at  approximately  200  ns,  and  then  sharply  increased  and  became 
unbounded  at  approximately  470  ns.  As  noted  by  FaBarbera  and  Zikry  [58],  this  unbounded 
thennal  increase  is  an  indication  of  hot  spot  formation,  and  this  is  also  the  region  where  there 
was  maximum  shear  strains  and  temperature  buildup  as  well  as  viscous  sliding  of  the  estane 
binder.  The  normal  stress,  shear  stress,  and  pressure  increased  shortly  before  50  ns.  Then, 
when  the  plastic  defonnation  initiated,  the  nonnal  stress,  shear  stress,  and  pressure  continued 
to  increase  significantly.  This  large  increase  in  temperatures  leads  to  thennal  decomposition, 
unbounded  temperature  buildup,  and  hot  spot  formation. 

The  experiments  of  Bruckman  and  Guillet  [48]  and  Field  et  al.  [59]  indicate  that  a  hot 
spot  can  fonn  or  initiate  when  the  heat  generated  due  to  decomposition  exceeds  the  heat 
dissipated  to  the  sunoundings,  which  results  in  run-away  or  unbounded  temperatures.  As  also 
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analyzed  by  LaBarbera  and  Zikry  [58],  the  heating  due  to  plastic  work  in  RDX  crystals  is 
substantial,  and  it  can  result  in  hot  spot  formation  for  dynamic  pressure  loadings,  and 
maximum  localized  thermal  accumulation  (Figure  5.4a).  The  mechanism  for  hot  spot 
initiation  for  the  amorphous  elastic  PCTFE  binder  is  heating  due  to  accumulated  plasticity  at 
the  peripheries  of  the  voids,  within  the  RDX  crystals,  where  the  shear  strain  localizes  [58]. 
However,  the  mechanism  for  hot  spot  initiation  in  RDX  aggregates  with  a  viscoelastic  estane 
binder  was  heating  due  to  plasticity  at  the  crystal  edges  associated  with  crystal  sliding  due  to 
the  viscoelastic  shearing  of  the  polymer  binder.  Grain  sliding  induced  regions  of  plastic 
deformation  at  the  edges  of  the  RDX  crystals,  which  resulted  in  hot  spot  formation  at  the 
edges  of  the  RDX  crystals.  Therefore,  binder  composition  significantly  influenced  the 
mechanism  for  hotspot  formation,  where  hot  spot  fonnation  occurs  by  shear  strain 
localization  at  void  peripheries  with  a  PCTFE  binder  or  occurred  by  shear  strain  localization 
at  crystal  edges  due  to  crystal  sliding  with  the  estane  binder.  Crystal  sliding  is  also  affected 
by  the  thermo-mechanical  incompatibility  with  RDX  crystals.  Dislocations  at  the  RDX- 
estane  interface  cannot  transmit  from  the  RDX  crystal  into  the  estane  binder  and  therefore 
dislocations  pileup  at  the  crystal  interface.  Pileups  can  result  in  an  avalanche  release  of  the 
dislocations  and  potentially  form  a  micro-crack  at  the  RDX-estane  interface  as  postulated  by 
Annstrong  and  Elban  [1].  These  dislocation-density  buildups  result  in  a  buildup  of  normal 
stress  and  flow  stress,  which  activates  slip  behavior,  and  then  results  in  increased  plastic 
accumulation  at  the  RDX-estane  interfaces. 
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5.2.3  Dynamic  Compressive  Pressure  Loading  of  1  GPa 

For  a  dynamic  compressive  loading  of  1  GPa,  at  a  time  of  .346  ps,  the  temperatures 

accumulated  at  the  crystal  edges  (Figure  5.9a)  due  to  grain  sliding.  The  maximum 
nonnalized  temperature  in  the  aggregate  was  1.97,  for  a  maximum  temperature  rise  of  283  K 
over  a  time  span  of  .346  ps,  which  resulted  in  hot  spot  formation  in  the  region  of  largest 
temperature  rise,  which  is  identified  in  Figure  5.9a.  The  nonnalized  normal  stress  (Figure 
5.9b)  was  generally  compressive  in  the  aggregate  with  a  minimum  normalized  value  of  -7. 
Specific  regions  of  tensile  normal  stress  occurred  both  in  the  polymer  binder  and  near  the 
edges  of  the  RDX  crystals,  with  a  maximum  nonnalized  value  of  1 .  This  is  the  result  of 
dynamic  wave  reflections  and  stress  incompatibilities  due  to  the  RDX-estane  interfaces.  The 
lateral  stress  (Figure  5.9c)  was  generally  compressive  with  a  maximum  compressive 
nonnalized  value  of  -4  that  occuned  within  different  RDX  crystals.  The  maximum 
nonnalized  value  of  tensile  lateral  stress  was  3,  and  occurred  at  the  interfaces  of  the  binder 
and  the  RDX  crystals.  The  maximum  values  of  shear  stress,  also  occuned  at  the  crystal- 
binder  interfaces,  and  the  maximum  value  was  -4  (Figure  5.9d).  The  maximum  normal  stress 
occurs  at  the  interfaces  due  to  the  viscoelastic  formulation  for  the  estane  binder  while  the 
lateral  stress  maximum  occurs  at  the  top  of  the  aggregate  where  the  crystals  are  under  tensile 
lateral  stress. 

The  energy  dissipated  (Figure  5.10a)  in  the  aggregate  was  due  to  plasticity  as  well  as 
viscous  effects  in  the  polymer  with  a  maximum  value  of  1 . 1E9  joules.  The  hydrostatic 
pressure  (Figure  5.10b)  had  a  maximum  tensile  nonnalized  value  of  1  and  a  maximum 
compressive  value  of  -3.5.  All  of  these  values  occurred  at  the  RDX-  binder  interfaces.  The 
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plastic,  inelastic,  shear  slip  (Figure  5.10c)  occurred  at  the  crystal  edges  and  had  a  maximum 
value  of  65%. 


The  mobile  dislocation  density  for  slip  system  one  (Figure  5.11a)  had  a  maximum 
nonnalized  value  of  3500.  The  mobile  dislocation  density  for  slip  system  two  (Figure  5.11b) 
had  a  maximum  nonnalized  value  of  260,  which  is  an  order  of  magnitude  less  than  slip 
system  one.  The  mobile  dislocation  density  for  slip  system  three  (Figure  5.11c)  had  a 
maximum  nonnalized  value  of  13,  which  is  an  order  of  magnitude  lower  than  slip  system 
two.  As  shown  in  LaBarbera  and  Zikry  [58],  low  random  GB  misorientations  cause 
substantial  slip  activity  to  occur  on  slip  system  one  with  less  substantial  activity  occurring  on 
slip  systems  two  and  three.  The  immobile  dislocation  density  for  slip  system  one  (Figure 
5.12a)  had  a  maximum  normalized  value  of  350,  while  slip  system  two  and  three  (Figure 
5. 12(b,c))  had  maximum  normalized  values  of  2.4  and  1 .2,  respectively. 

The  temperature,  nonnal  stress,  shear  stress,  and  pressure  evolution,  as  a  function  of 
time  history,  for  the  element  that  had  the  highest  temperature  rise  and  subsequently  lead  to 
hot  spot  fonnation,  is  shown  in  Figure  5.13.  This  element  was  located  on  the  edge  of  the 
specified  crystal  as  shown  in  Figure  5.9a.  The  temperature  was  constant  until  plastic 
deformation  occurred  at  approximately  210  ns  after  application  of  the  load,  and  as  a  result  the 
temperature  increased.  At  approximately  290  ns,  the  temperature  increased  sharply  and  within 
60  ns  the  temperature  became  unbounded  due  to  the  dominant  effects  of  thermal 
decomposition.  The  normal  stress  becomes  highly  compressive  at  approximately  50  ns  after 
the  initial  application  of  loading  and  remains  highly  compressive  until  the  temperature 
increased  sharply,  at  which  time  it  began  to  unload.  The  shear  stress  had  a  relatively  low 
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value  of  approximately  - 1 ,  until  plastic  defonnation  occurred  and  the  shear  stress  sharply 
increased  to  approximately  -4.5  at  a  time  of  340  ns.  The  hydrostatic  pressure  was  initially 
compressive  until  the  temperature  became  unbounded,  at  which  time  the  pressure 
transitioned  from  compressive  to  tensile  in  the  element  where  the  hot  spot  resulted.  This 
indicates  that  as  hot  spot  fonnation  occurs  there  is  a  complex  interaction  with  the 
microstructure  that  results  in  an  unbounded  temperature,  an  unloading  of  normal  stress,  a 
buildup  of  shear  stress,  and  a  transition  of  pressure  from  compressive  to  tensile.  This  also 
indicates  that  the  region  where  the  hot  spot  fonnation  occurs  is  dominated  by  a  buildup  of 
shear  stress. 

As  noted  earlier  for  the  1  GPa  tensile  pressure  loading  case,  the  primary  mechanism 
for  hot  spot  formation  was  inelastic  defonnation  from  the  resulting  dislocation  motion  due  to 
energy  dissipation  associated  with  crystal  sliding.  Crystal  sliding  was  facilitated  by  the 
application  of  a  high  dynamic  pressure  loading,  which  resulted  in  shear  deformation  of  the 
polymer  binder,  and  therefore  caused  crystal  sliding  and  crystal-  estane  interaction.  For  the 
tensile  loading  case,  the  time  to  hotspot  fonnation  was  greater  in  comparison  with  the 
compressive  loading  case,  which  is  due  to  the  incompressibility  of  the  polymer  binder,  in  the 
compression  case,  resulting  in  RDX  plastic  deformation  at  an  earlier  time. 

The  deformation  that  occuned  as  a  result  of  the  viscoelastic  polymer  binder  led  to 
complex  interactions  at  RDX  crystal  to  binder  interfaces.  Figure  5.14b  shows  the  spatial 
temperature  and  pressure  distribution  across  the  RDX  crystal  -  binder  interface  that  is 
denoted  in  Figure  5.14a  by  the  thick,  dark  straight  line.  The  temperature  had  accumulated  at 
the  edges  of  the  RDX  crystal  that  is  in  contact  with  the  polymer  binder.  As  a  result,  the 
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pressure  is  highly  compressive  in  the  RDX  crystals,  and  it  increases  until  the  polymer  binder 
is  encountered.  This  increase  in  pressure  is  due  to  the  incompressibility  of  the  polymer 
material.  This  is  the  main  mechanism  that  results  in  crystal  sliding  due  to  the  ease  of  polymer 
shearing  and  incompressibility  of  the  polymer.  This  indicates  that  relatively  small  regions  of 
the  polymer,  which  are  constrained  between  two  neighboring  crystals,  are  unable  to  compress 
but  are  easily  sheared,  and  this  facilitated  RDX  crystal  sliding.  This  also  indicates  that  as  the 
crystal  slides  due  to  incompressibility,  the  large  pressures  resulted  in  plastic  defonnation  at 
the  crystal  edges,  which  led  to  temperature  accumulation,  and  this  subsequently  can  lead  to 
hot  spot  fonnation.  For  the  tensile  loading  case,  the  regions  of  hydrostatic  compression 
occurred  in  the  vertical  regions  of  the  binder  ligaments  (Figure  5.5b)  and  for  compressive 
loading  the  hydrostatic  compression  occurred  in  the  horizontal  binder  regions  (Figure  5. 10b). 
For  both  the  compressive  and  tensile  loading  cases,  the  hot  spot  fonnation  mechanism  is  due 
to  hydrostatic  pressures,  which  resulted  in  plastic  defonnation  and  unbounded  temperatures. 
Due  to  the  tensile  pressure,  there  is  a  material  competition  between  plastic  deformation  and 
debonding  at  the  RDX-estane  interface  that  could  lead  to  further  incompatibilities  at  binder- 
RDX  interfaces. 

5.3  Summary 

A  dislocation-density  based  crystalline  plasticity,  viscoelasticity,  and  specialized 
finite-element  fonnulations  were  used  to  investigate  the  mechanisms  for  hot  spot  fonnation 
due  to  dynamic  thenno-mechanical  loading  conditions  in  RDX-estane  energetic  aggregates. 
For  a  1  GPa  dynamic  tensile  pressure  loading,  the  temperature  accumulated  at  the  periphery 
of  a  specific  crystal,  where  the  shear  strains  had  localized  at  the  estane-RDX  interfaces.  The 
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temperature  accumulation  resulted  in  hot  spot  formation  at  the  edges  of  crystals.  The  energy 
due  to  plastic  defonnation  localized  at  the  peripheries  of  the  crystals  and  viscous  dissipation 
was  localized  in  the  binder  regions,  which  had  extensively  sheared,  and  resulted  in  crystal 
sliding.  Plastic  deformation  was  dominated  by  one  dominant  slip  system  due  to  the  random 
low  GB  misorientations. 

For  the  1  GPa  dynamic  compressive  loading,  the  temperature  localized  at  the 
peripheries  of  the  crystals,  where  extensive  plastic  defonnation  had  accumulated.  Plastic 
deformation  was  facilitated  by  crystal  sliding  due  to  extensive  shearing  and  the 
incompressibility  of  the  estane  binder.  At  RDX-estane  interfaces  that  had  extensively 
sheared,  the  temperature  accumulated  and  the  magnitude  of  temperature  and  pressure  buildup 
significantly  increased  as  the  interface  was  approached.  This  increase  in  hydrostatic  pressure 
was  due  to  a  combination  of  the  relative  ease  of  polymer  shearing  and  incompressibility. 
These  regions  of  incompressibility  in  the  binder  facilitated  crystalline  deformations  that 
resulted  in  localized  plastic  deformation  and  hot  spot  formation. 

The  mechanisms  for  hot  spot  fonnation,  for  both  compressive  and  tensile  loading 
conditions,  for  a  plastically  bonded  explosive  at  an  initial  temperature  that  is  significantly 
above  the  glass  transition  temperature,  is  RDX  inelastic  defonnations  at  the  crystal  edges  and 
is  facilitated  by  crystal  sliding  at  the  RDX-estane  interfaces.  This  results  in  localized  plastic 
deformation  at  the  periphery  of  the  crystals.  This  is  in  contrast  to  a  plastically  bonded 
explosive  at  a  temperature  significantly  less  than  the  glass  transition  temperature  (see,  for 
example  [58]),  where  shear  strain  localization  at  the  peripheries  of  the  voids  results  in  hot 
spot  formation.  Incompatibilities  between  the  RDX  crystal  and  the  estane  binder  exist  due  to 
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dislocation-density  accumulation  at  crystalline  interfaces,  and  this  results  in  a  buildup  of 
inelastic  defonnation  at  the  RDX-estane  interface. 

Results  from  this  study  indicate  that  the  interrelated  effects  of  dynamic  loading, 
polymer  binder  response,  and  void  porosity  significantly  affect  hot  spot  formation 
mechanisms.  Unbounded  temperatures  result  from  localized  plastic  defonnation  and  heat 
released  due  to  thennal  decomposition,  which  was  thermally  activated  by  the  heat  generated 
due  to  plastic  deformation  and  viscoelastic  sliding.  Different  polymer  binders  with  different 
glass  transition  temperatures  will  affect  the  thenno-mechanisms  for  hot  spot  formation  due  to 
thenno-mechanical  interfacial  incompatibilities  with  RDX  crystals. 
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5.4  Tables  and  Figures 


Table  5.1:  Material  Properties  of  RDX  and  Estane  binder 


RDX 

Estane  Binder 

Density  (!  ) 

1 820 

TT9B 

1  ! 

Yield  Stress  (MPa) 

580 

— 

Elastic  Modulus  (GPa) 

18.4 

— 

Poisson  Ratio 

0.22 

.499 

f 

Thermal  Conductivity  (  '  ) 

0.29 

0.14 

f 

Specific  Heat  (cp)  (— ) 

!  "  f 

1260 

1500 

Rate  Sensitivity  Parameter 

.01 

— 

Initial  Mobile  Dislocation 

Density  (m"  ) 

lelO 

— 

Initial  Immobile  Dislocation 
Density  (m'2) 

lel2 

— 
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Figure  5.1 :  The  RDX-Estane  aggregate  with  loading,  boundary  conditions,  and  void 

distribution. 
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Figure  5.2:  The  Estane  model  with  loading  and  boundary  conditions  used  for  model 

validation. 
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Figure  5.3:  (a)  Quasi-static  strain  rate  of  0.005  s'1,  (b)  dynamic  strain  rate  of  2000  s'1 
comparison  of  viscoelastic  model  with  experimental  data  for  estane. 
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Figure  5.4:  (a)  Temperature,  (b)  normal  stress,  (c)  lateral  stress,  and  (d)  shear  stress  for  a  1 
GPA  tensile  loading  with  low  angle  GB  misorientations.  (Spatial  dimensions  are  in  units  of 

meters) 
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Figure  5.5:  (a)  Dissipated  energy,  (b)  pressure,  and  (c)  plastic  shear  slip  for  a  1  GPA  tensile 
loading  with  low  angle  GB  misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  5.6:  (a)  Mobile  dislocation  density  l,  (b)  mobile  dislocation  density  2,  and  (c)  mobile 
dislocation  density  3  for  a  1  GPA  tensile  loading  with  low  angle  GB  misorientations.  (Spatial 

dimensions  are  in  units  of  meters) 


69 


(a) 


(b) 


(010)  |001) 


00012 

0001 

oocce 

00X6 

00004 

OOCOP 

0 


Immobile  Dislocation 
Density  1 

800 
700 
600 
500 
400 
300 
200 
100 


Immobile  Dislocation 
Density  3 

!6 

5.5 
5 

4.5 
4 

3.5 
3 

2.5 
2 

1.5 


00OCP  00001  OCOCC  O0CC8  0001  00012  00014 


0001? 


I  ...  I 

'o  ooocp  oocw  oocco  oooce  oooi  oooi?  00014 

X 


Immobile  Dislocation 
Density  2 

100 
90 
80 
70 
60 
50 
40 
X 
20 
10 


(021) [ioo| 


0001? 


oat* 


00001 


(021)1100) 


Figure  5.7:  (a)  Immobile  dislocation  density  1,  (b)  immobile  dislocation  density  2,  and  (c) 
immobile  dislocation  density  3  for  a  1  GPA  tensile  loading  with  low  angle  random  GB 
misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  5.8:  Temperature,  normal  stress,  shear  stress,  and  pressure  time  history  for  the  region 
where  a  hotspot  fonnation  occurred  for  a  1  GPA  tensile  loading  with  low  angle  GB 

misorientations. 
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Stress  and  Pressure  (Normalized  by  580  MPa) 


Figure  5.9:  (a)  Temperature,  (b)  normal  stress,  (c)  lateral  stress,  and  (d)  shear  stress  for  a  1 
GPA  compressive  loading  with  low  angle  GB  misorientations.  (Spatial  dimensions  are  in 

units  of  meters) 
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Figure  5.10:  (a)  Dissipated  energy,  (b)  pressure,  and  (c)  plastic  shear  slip  for  a  1  GPA 


compressive  loading  with  low  angle  GB  misorientations.  (Spatial  dimensions  are  in  units  of 

meters) 
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Mobile  Dislocation 


Figure  5.11:  (a)  Mobile  dislocation  density  1,  (b)  mobile  dislocation  density  2,  and  (c)  mobile 
dislocation  density  3  for  a  1  GPA  compressive  loading  with  low  angle  GB  misorientations. 

(Spatial  dimensions  are  in  units  of  meters) 
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Figure  5.12:  (a)  Immobile  dislocation  density  1,  (b)  immobile  dislocation  density  2,  and  (c) 
immobile  dislocation  density  3  for  a  1  GPA  compressive  loading  with  low  angle  GB 
misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  5.13:  Temperature,  normal  stress,  shear  stress,  and  pressure  time  history  for  the  region 
where  a  hotspot  formation  occurred  for  a  1  GPA  compressive  loading  with  low  angle  GB 

misorientations. 
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Figure  5.14:  (a)  The  spatial  location  where  temperature  and  pressure  data  were  extracted 
(specifically  along  the  straight  line  within  the  circled  region)  and  (b)  is  the  spatial  distribution 
for  temperature  and  pressure  across  the  RDX-Estane  interface. 
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Pressure  (Normalized  by  580  MPa) 


CHAPTER  6:  Dynamic  Fracture  and  Local  Failure  Mechanisms  in  Heterogeneous 
RDX-Estane  Energetic  Aggregates 

6.1  Introduction 

The  heterogeneous  crystalline-amorphous  structure  of  energetic  materials,  the  defects 
that  span  different  spatial  scales,  the  fracture  susceptibility  of  these  crystalline-amorphous 
aggregates,  and  the  dynamic  thermo-mechanical  loading  conditions  underscore  that  a 
fundamental  understanding  and  prediction  of  local  failure  initiation  behavior  is  needed.  The 
objective  of  this  investigation,  therefore,  is  to  formulate  a  predictive  framework  that  accounts 
for  the  interrelated  micro  structural  effects  of  crystalline  plasticity,  polymer  binder  relaxation, 
crystalline-amorphous  fracture,  and  inherent  microstructural  defects  and  features,  such  as 
void  distribution,  dislocation-density  evolution,  and  GBs  on  dynamic  crack  nucleation  and 
growth  in  aggregates  subjected  to  high  strain  rate  loads.  A  recently  developed  dislocation- 
density  based  crystalline  plasticity  [35]  approach,  which  has  been  coupled  to  a  finite 
viscoelasticity  formulation,  dynamic  fracture  formulation,  and  dynamic  finite-element 
techniques  will  be  used  to  understand  local  failure  mechanisms  in  RDX-estane  aggregates. 

6.2  Results 

Computational  analyses  were  undertaken  to  investigate  and  predict  the  effects  of 
dynamic  strain  rate  loading  conditions  on  crack  nucleation  behavior  and  local  failure  for  a 
RDX  crystalline  aggregate  with  a  viscoelastic  polymer  binder.  The  finite-element  aggregate 
of  100  grains  of  RDX  crystals  with  an  average  grain  size  of  200  pm  for  an  overall  aggregate 
size  of  2mm  x  2mm  (Figure  6.1)  consisted  of  13,382  quadrilateral  elements. 
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The  RDX  crystalline  morphology  was  generated  utilizing  a  Voronoi  tessellation  scheme. 
Based  on  Gallagher  et  al.  [44],  RDX  grains  were  assumed  to  have  the  following  three 
independent  slip  systems  -  (010)  [001],  (021)  [100],  (02-1)  [100]  -  with  material  properties 
given  by  Armstrong  and  Elban[l]  and  Annapragada  et  al.  [50]  in  Table  6.1.  A  polymer 
binder  of  estane,  with  a  volume  fraction  of  approximately  19%,  was  randomly  distributed 
between  the  RDX  grains.  Material  properties  for  the  estane  polymer  binder  were  based  on 
measurements  by  Mas  et  al.  [39]  and  Barua  and  Zhou  [23],  as  shown  in  Table  6. 1 . 

Porosity  in  RDX  crystals  can  be  controlled  by  varying  crystal  processing  parameters, 
which  can  result  in  a  vastly  different  combination  of  void  morphologies,  sizes,  and 
distributions  [7,  46],  In  this  investigation,  circular  voids,  with  an  average  diameter  of  15  pm 
and  an  overall  porosity  volume  fraction  of  2.5%  were  used  (Figure  6. 1),  which  is  consistent 
with  energetic  material  microstructures  [7,  46,  49,  51,  52].  The  domain  of  influence  radius 
for  fracture  nucleation  was  set  to  15  pm.  The  fracture  stress  was  set  to  approximately  three 
times  the  yield  stress  of  RDX  for  both  RDX  and  estane  binder.  Fixed  boundary  conditions 
were  applied  to  the  model  bottom  for  two-dimensional  plane-strain  loading  conditions.  A 
tensile  velocity  was  applied  on  the  top  surface,  as  shown  in  Figure  6.1,  for  dynamic  nominal 
strain  rates  of  1,000  s'1  and  10,000  s'1,  with  random  GB  misorientations  of  less  than  15°. 

6.2.1 1,000  s'1  Nominal  Strain  Rate 

For  a  nominal  strain  rate  of  1,000  s'1  at  a  nominal  strain  of  9%,  the  maximum 
nonnalized  normal  stress  (Figure  6.2a)  was  5  (stresses  were  nonnalized  by  the  static  yield 
strength  of  RDX),  which  occurred  at  the  interface  of  the  RDX-estane  where  a  crack  had 
nucleated  from  the  adjacent  void  periphery  and  subsequently  propagated  towards  the 
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interface  and  then  was  arrested.  The  nonnalized  nonnal  stress  was  a  minimum  of  -2,  which 
occurred  in  the  thin  ligaments  of  binder  that  were  under  hydrostatic  compression.  From  the 
normal  stress  spatial  distribution  in  Figure  6.2a,  multiple  cracks  have  nucleated  at  the 
peripheries  of  the  voids  where  the  cleavage  stresses  were  large  and  subsequently  several  of 
these  cracks  propagated  to  the  estane  binder  where  the  cracks  were  arrested,  which  is 
consistent  with  experimental  results  [29-31].  The  lateral  stress  (Figure  6.2b)  had  a  maximum 
nonnalized  value  of  2.5,  which  occurred  in  the  horizontal  ligaments  of  the  binder  due  to  large 
tensile  pressures  in  the  binder  and  had  a  minimum  of  -4.  The  shear  stress  (Figure  6.2c)  had  a 
maximum  nonnalized  value  of  1 .2  and  a  minimum  of -1.8,  which  both  occurred  at  the  corners 
of  specific  crystals  that  had  undergone  extensive  shearing  (Figure  3  c)  in  the  polymer  binder. 
The  hydrostatic  pressure  (Figure  6.2d)  had  a  maximum  normalized  value  of  2.2  and  a 
minimum  of  -1.6  which  both  occurred  in  the  binder  ligaments  due  to  the  bulk 
incompressibility  of  the  binder  under  hydrostatic  conditions. 

The  temperature  spatial  distribution  (Figure  6.3a)  had  an  accumulation  of 
temperatures  at  the  peripheries  of  the  crystals  for  a  maximum  nonnalized  value  of  1 .22 
(temperature  was  nonnalized  by  the  initial  temperature  of  293  K).  This  buildup  of 
temperature  is  due  to  the  combined  effects  plastic  work  dissipation  and  thermal 
decomposition.  The  maximum  buildup  of  thennal  accumulation  was  located  at  the  periphery 
of  the  specific  grain  shown  in  Figure  6.3c  and  had  localized  due  to  the  competition  between 
crack  tip  propagation  and  crack  arrest  at  the  interface  which  resulted  in  a  large  buildup  of 
temperature.  The  interaction  of  crack-tip  propagation  and  crack  arrest  is  consistent  with 
experimental  results[29-3 1],  The  shear  slip,  which  is  the  sum  of  the  plastic  strains  on  all  slip 
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systems,  had  localized  in  similar  regions  as  the  temperature  accumulations  had  and  was  a 
maximum  of  22%.  Figure  6.3d  shows  the  dissipated  energy  due  to  plastic  work  and  viscous 
dissipation,  where  the  maximum  energy  dissipated  was  approximately  2E8  joules  and 
occurred  at  the  same  spatial  location  as  the  largest  temperature  accumulation. 

The  mobile  dislocation  densities  for  slip  system  one  (Figure  6.4a),  (010)  [001],  which 
was  the  most  active  slip  system,  had  a  maximum  nonnalized  accumulation  of  1000.  The 
dislocation  densities  were  nonnalized  by  the  initial  density  values.  The  buildups  occurred  at 
the  peripheries  of  the  crystals  where  the  shear  strains  had  localized.  The  immobile  dislocation 
densities  for  slip  system  one  (Figure  6.4b),  (010)  [001],  had  a  maximum  normalized  value  of 
100  which  also  occuned  at  the  periphery  of  a  specific  crystal. 

Due  to  the  material  competition  between  crack  propagation  and  crack  arrest  at  the 
interface,  a  large  buildup  of  temperature  occurred  at  the  interface  as  shown  in  Figure  6.3c. 
Therefore,  the  nonnal  stress,  pressure,  shear  stress,  lateral  stress,  mobile  dislocation  density, 
and  immobile  dislocation  density  time  history  evolution  for  the  interface-crack  interaction 
region  shown  in  Figure  6.3c  is  given  in  Figure  6.5.  The  nonnal  stress,  pressure,  shear  stress, 
and  lateral  stress  increased  linearly,  in  magnitude,  for  a  time  of  approximately  85  pscc.  At 
which  time  a  crack  nucleated  at  the  nearby  void  periphery  and  subsequently  began  to 
propagate  towards  the  interface  at  a  time  of  approximately  90  usee,  which  marked  a  sharp 
increase  in  the  nonnal  stress,  pressure,  shear  stress,  and  lateral  stress.  The  mobile  and 
immobile  dislocation  densities  were  constant  until  a  time  of  approximately  90  pscc,  when  the 
crack  propagated  to  the  nearby  interface  and  then  the  mobile  dislocation  densities  increased 
sharply  to  260  (the  dislocation  densities  were  normalized  by  the  initial  density  values)  and 
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the  immobile  dislocation  densities  increased  sharply  to  approximately  25.  The  sharp  increase 
in  normal  stress,  pressure,  shear  stress,  lateral  stress,  mobile  dislocation  density,  and 
immobile  dislocation  density,  at  a  time  of  approximately  90  usee  is  due  to  the  interaction  of 
the  crack  tip  propagating  to  the  interface  and  impeding  further  propagation  due  to  the  low 
stress  barrier  in  the  polymer  binder,  this  subsequently  resulted  in  a  buildup  of  mobile  and 
immobile  dislocation  densities  at  the  interface.  This  crack  tip  interaction  with  interfacial 
regions  underscores  how  crack  arrest  can  be  affected  by  the  low  stresses  in  the  polymer 
regions  as  experimentally  observed  [27-31], 

The  stress-strain  curve  for  a  nominal  strain  rate  of  1,000  s'1  is  shown  in  Figure  6.6. 
The  stress  linearly  increased  with  strain  until  a  nominal  strain  of  7.5%,  at  which  point  the 
stress  is  approximately  0.95  (stress  has  been  nonnalized  by  the  static  yield  stress  of  RDX) 
and  begins  to  unload  due  to  localized  plasticity  and  crack  nucleation,  which  are  both  local 
failure  mechanisms  within  the  aggregate.  Small  dynamic  oscillations  in  the  stress-strain 
curve  occurred  at  low  nominal  strains,  which  dampened  out  after  approximately  at  a  nominal 
strain  of  1%.  The  aggregate  had  a  maximum  stress  of  1  at  a  nominal  strain  of  8.5%.  At  a 
nominal  strain  of  7.5%,  multiple  cracks  nucleated  and  propagated,  and  this  resulted  in  an 
unstable  failure  at  a  nominal  strain  of  9%.  Therefore,  the  competing  effects  of  plasticity, 
thennal  decomposition,  crack  nucleation,  dislocations,  and  voids  led  to  unstable  failure,  as 
indicated  by  the  stress  unloading  (Figure  6.5)  where  the  local  stresses  and  dislocation 
densities  sharply  increased  as  the  crack  impinged  on  the  RDX-estane  interface. 
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6.2.2  10,000  s' 1  Nominal  Strain  Rate 

For  a  nominal  strain  rate  of  10,000  s'1  at  a  nominal  strain  of  9.4%  the  nonnal  stress 
(Figure  6.7a)  had  a  maximum  normalized  value  of  5.5  that  occurred  at  the  RDX  crystal- 
estane  interface.  Crack  nucleated  from  voids,  and  they  subsequently  propagated  towards  the 
estane  interface.  This  is  consistent  with  experimental  observations  that  cracks  can  be 
impeded  by  viscoelastic  or  amorphous  regions  [29-31],  The  buildup  of  nonnal  stress  at  the 
interface  is  also  similar  with  the  1,000  s'1  nominal  strain  rate  case.  The  only  difference  is  that 
the  normal  stresses  were  higher  for  the  10,000  s'1  case,  which  is  due  to  dynamic  strain  rate 
hardening.  The  nonnal  stress  had  a  minimum  nonnalized  value  of -1.5,  which  occurred  in  the 
vertical  ligaments  of  the  binder  that  were  under  hydrostatic  compression  (Figure  6.7d).  The 
lateral  stress  (Figure  6.7b)  had  a  maximum  nonnalized  value  of  2.5,  which  occurred  in  the 
horizontal  ligaments  of  the  binder  that  were  generally  under  hydrostatic  tensile  pressure  as 
evidenced  from  Figure  6.7d.  The  lateral  stress  had  a  minimum  nonnalized  value  of  -5.5.  The 
shear  stress  (Figure  6.7c)  had  a  maximum  normalized  value  of  1.2  and  a  minimum 
nonnalized  value  of  -1.4,  and  these  maximum  values  occurred  at  the  peripheral  regions  of  the 
RDX  crystals  where  the  adjacent  estane  binder  had  experience  extensive  shearing.  The 
hydrostatic  pressure  (Figure  6.7d)  had  a  maximum  normalized  value  of  2.5,  which  occurred 
in  the  horizontal  ligaments  of  the  binder  due  to  the  large  bulk  modulus  associated  with  the 
polymer  binder.  The  hydrostatic  pressure  had  a  minimum  nonnalized  value  -1.5  that  occurred 
in  the  vertical  ligaments  of  the  binder  that  were  under  hydrostatic  compression  due  to  crystal 
sliding. 
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The  temperature  (Figure  6.8a)  accumulated  at  the  peripheries  of  the  voids  and  crystals 
as  well  as  in  regions  where  crack-tip  interaction  with  the  interface  had  occurred  with  a 
maximum  normalized  value  of  1 .45  occurring  at  the  crack-tip  interaction  region  shown  in 
(Figure  6.8c).  The  plastic  shear  slip  (Figure  6.8b)  accumulated  a  maximum  value  of  45%, 
which  occurred  at  the  crack-tip  interaction  region  shown  in  (Figure  6.8c).  Figure  6.8c  shows 
a  zoomed  in  view  of  the  crack-tip  interaction  with  the  polymer  binder  and  the  resultant 
temperature  rise.  The  dissipated  energy  (Figure  6.8d)  had  a  maximum  value  of  4.5E8  joules 
of  energy  dissipated,  which  occurred  at  the  crack-tip  interaction  region.  Due  to  the 
localization  of  temperature,  shear  slip,  and  dissipated  energy  at  the  crack-tip  interaction 
region,  the  system  is  dominated  by  crack  nucleation  at  the  void  peripheries  and  the 
subsequent  propagation  to  the  binder  where  hotspot  formation  occurred  as  detailed  by 
LaBarbera  and  Zikry  [58], 

The  mobile  dislocation  densities  for  slip  system  one  (Figure  6.9a),  (010)  [001],  which 
was  the  most  active  slip  system,  had  a  maximum  nonnalized  value  of  4200.  The  immobile 
dislocation  densities  for  slip  system  one  (Figure  6.9a),  (010)  [001],  had  a  maximum 
nonnalized  accumulation  of 420.  The  maximum  nonnalized  values  for  the  mobile  and 
immobile  dislocations  densities  occuned  at  the  crack-tip  interaction  region  as  shown  in 
Figure  6.8c. 

Figure  6.10  shows  the  temporal  evolution  of  the  nonnal  stress,  lateral  stress,  shear 
stress,  hydrostatic  pressure,  mobile  dislocation  density,  and  immobile  dislocation  density 
from  the  critical  interface  region  as  shown  in  Figure  6.8c.  The  nonnal  stress,  lateral  stress, 
shear  stress,  and  hydrostatic  pressure  all  increased,  in  magnitude,  up  to  a  nominal  strain  of 
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8.5%.  During  this  increase,  there  were  wave  reflections.  At  a  nominal  strain  of  8.5%,  the 
crack  that  has  nucleated  from  the  adjacent  void  (Figure  6.8c)  has  propagated  into  the  critical 
fracture  region  and  results  in  a  sharp  increase  in  the  nonnal  stress,  lateral  stress,  shear  stress, 
and  hydrostatic  pressure.  During  this  sharp  increase  in  normal  stress,  lateral  stress,  shear 
stress,  and  hydrostatic  pressure,  the  mobile  and  immobile  dislocation  densities  sharply 
increase.  At  a  nominal  strain  of  approximately  9.2%,  local  failure  occurred  and  subsequent 
unloading  of  stresses  began  as  indicated  by  the  sharp  drop  in  normal  stress,  lateral  stress, 
shear  stress,  and  hydrostatic  pressure.  Therefore,  the  crack  propagated  into  the  region  and 
encountered  the  interface.  Then  subsequently,  increased  plasticity  at  the  interface  resulted  in 
unstable  local  failure  of  the  aggregate  as  indicated  by  the  unloaded  stresses  (Figure  6.10). 

The  stress-strain  curve  for  the  overall  aggregate  response  is  shown  in  Figure  6. 1 1 .  At 
a  nominal  strain  of  approximately  1.5%,  initial  high  frequency  dynamic  wave  reflections 
dampened  out.  Up  to  9%  nominal  strain,  the  stress  increased  with  strain  while  also  exhibiting 
dynamic  wave  reflections.  At  a  nominal  strain  of  9%,  local  failure  and  dynamic  localized 
plasticity  began  to  occur  which  resulted  in  a  subsequent  lowering  in  the  stress  as  the  strain 
increased.  The  maximum  global  stress  was  approximately  1.1  at  a  nominal  strain  of  9%  , 
which  was  higher  than  the  maximum  normalized  stress  of  1  for  the  1,000  s'1  nominal  strain 
rate  case,  and  this  is  due  to  dynamic  strain-rate  hardening  . 

For  both  the  1,000  s'1  and  10,000  s'1  nominal  strain  rate  cases,  crack  nucleation 
occurred  primarily  at  the  peripheries  of  the  voids  where  the  stresses  on  the  cleavage  planes 
exceeded  the  local  fracture  stress.  Once  crack  nucleation  occurred,  the  cracks  propagated 
towards  the  binder  where  the  cracks  were  arrested.  This  is  due  to  the  propagating  crack 
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interacting  with  the  interrelated  effects  of  plasticity,  temperature,  dislocation  densities,  and 
shear  slip.  Crack  bridging  can  also  occur  across  the  viscoelastic  estane  binder  as  noted  by 
several  experimental  investigations  [29-3 1  ] . 

6.3  Summary 

The  effects  of  local  failure  initiation  through  crack  nucleation  and  propagation  were 
investigated  for  nominal  strain  rates  of  1,000  s'1  and  10,000  s'1.  Cracks  nucleated  at  void 
peripheries  and  propagated  to  the  estane  binder.  Crack  nucleation  and  propagation  was 
affected  by  the  heterogeneous  nature  of  RDX  energetic  aggregate  and  the  inherent  interfaces 
within  the  material.  The  RDX-estane  interfaces  affected  crack  propagation  and  resulted  in 
crack  arrest  at  the  binder.  Thennal  buildups  were  affected  by  crack  nucleation  and 
propagation  and  the  RDX-estane  interfaces.  Large  thermal  buildups  occurred  at  the  interfaces 
where  there  was  a  competition  between  crack  propagation  through  the  interface  and  crack 
arrest  at  the  interface. 

For  a  nominal  strain  rate  of  1,000  s'1,  the  nonnal  stresses  were  a  maximum  at  the 
crack-tip  interface  region.  The  shear  stress  was  a  maximum  at  the  RDX  crystal  peripheries 
where  extensive  shearing  had  occurred  in  the  adjacent  viscoelastic  estane  binder.  The  lateral 
stress  and  hydrostatic  tensile  pressure  were  a  maximum  in  the  horizontal  ligaments  of  the 
binder,  and  the  hydrostatic  compressive  pressure  was  a  maximum  in  vertical  ligaments  of  the 
binder  where  crystal  to  crystal  interaction  and  sliding  had  occurred.  The  temperature,  shear 
slip,  and  dissipated  energy  were  all  maximums  at  the  critical  interface  region  where  crack-tip 
interaction  occurred.  The  mobile  and  immobile  dislocation  densities  accumulated,  at  the 
peripheries  of  the  voids  and  crystals,  and  at  the  crack-tip  interaction  region.  The  temperature 
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buildup  was  dominated  by  the  competition  between  crack  propagation  through  the  RDX- 
estane  interface  and  crack  arrest  at  the  interface.  This  resulted  in  an  incompatibility,  for  crack 
propagation  through  the  interface,  which  affected  crack  growth.  For  a  nominal  strain  rate  of 
10,000  s'1,  the  results  were  similar  when  compared  with  a  strain  rate  of  1,000  s'1.  However, 
more  pronounced  dynamic  wave  reflections  were  shown  in  the  stress-strain  curve  for  a 
nominal  strain  rate  of  10,000  s'1  in  comparison  to  a  nominal  strain  rate  of  1,000  s'1. 

Results  from  this  study  indicate  that  the  interrelated  effects  of  dynamic  strain  rate 
loading,  polymer  binder  response,  void  porosity,  and  crack  nucleation/propagation 
significantly  affect  local  failure  mechanisms.  Voids  facilitated  the  nucleation  of  cracks, 
which  then  resulted  in  local  failure  of  the  aggregate  as  the  cracks  propagate  toward  the 
binder.  The  polymer  binder  affected  interfacial  incompatibilities  with  RDX  crystal  fracture 
and  resulted  in  local  failure  of  the  aggregate  due  to  a  buildup  of  inelastic  deformation  at  the 
interface  and  subsequent  failure. 
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6.4  Tables  and  Figures 


Table  6.1:  RDX  and  Estane  binder  Material  Properties 


RDX 

Estane  Binder 

Density  (!  ) 

1 820 

TT90 

I  ! 

Yield  Stress  (MPa) 

580 

— 

Elastic  Modulus  (GPa) 

18.4 

— 

Poisson  Ratio 

0.22 

.499 

V 

Thermal  Conductivity  (  '  ) 

0.29 

0.14 

1 

Specific  Heat  (cp)  (— -) 

!  "  ! 

1260 

1500 

Rate  Sensitivity  Parameter 

.01 

— 

Fracture  Stress  (GPa) 

1.8 

1.8 

Initial  Mobile  Dislocation 

Density  (m"  ) 

lelO 

— 

Initial  Immobile  Dislocation 
Density  (m"  ) 

lel2 

— 

88 


v(t) 


Figure  6. 1 :  The  RDX-estane  aggregate  with  loading  and  boundary  conditions,  void 


distribution,  and  grain  morphology. 
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Figure  6.2:  (a)  Normal  Stress,  (b)  lateral  stress,  (c)  shear  stress,  and  (d)  pressure  for  a 
nominal  strain  rate  loading  of  1,000  s'1  and  at  a  nominal  strain  of  9.0%.  (Spatial  dimensions 

are  in  units  of  meters) 
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Figure  6.3:  (a)  Temperature,  (b)  shear  slip,  (c)  zoomed  in  view  of  (a),  and  (d)  dissipated 
energy  for  a  nominal  strain  rate  loading  of  1,000  s'1  and  at  a  nominal  strain  of  9.0%.  (Spatial 

dimensions  are  in  units  of  meters) 
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Figure  6.4:  (a)  Mobile  and  (b)  immobile  dislocation  densities  for  the  most  active  slip  system 
for  a  nominal  strain  rate  loading  of  1,000  s'1  and  at  a  nominal  strain  of  9.0%.  (Spatial 

dimensions  are  in  units  of  meters) 
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Figure  6.5:  Normal  stress,  pressure,  shear  stress,  lateral  stress,  mobile  and  immobile 
dislocation  density,  for  slip  system  one,  time  history  for  the  interface-crack  interaction  region 

for  a  nominal  strain  rate  loading  of  1,000  s’1. 
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Figure  6.6:  Stress-Strain  curve  for  a  nominal  strain  rate  loading  of  1,000  s'1. 
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Figure  6.7:  (a)  Normal  Stress,  (b)  lateral  stress,  (c)  shear  stress,  and  (d)  pressure  for  a  nominal 
strain  rate  loading  of  10,000  s'1  and  at  a  nominal  strain  of  9.4%.  (Spatial  dimensions 

are  in  units  of  meters) 
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Figure  6.8:  (a)  Temperature,  (b)  shear  slip,  (c)  zoomed  in  view  of  (a),  and  (d)  dissipated 
energy  for  a  nominal  strain  rate  loading  of  10,000  s'1  and  at  a  nominal  strain  of  9.4%. 
(Spatial  dimensions  are  in  units  of  meters) 
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Figure  6.9:  (a)  Mobile  and  (b)  immobile  dislocation  densities  for  the  most  active  slip  system 
for  a  nominal  strain  rate  loading  of  10,000  s'1  and  at  a  nominal  strain  of  9.4%.  (Spatial 

dimensions  are  in  units  of  meters) 
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Figure  6.10:  Normal  stress,  pressure,  shear  stress,  lateral  stress,  mobile  and  immobile 
dislocation  density  time  history  for  the  interface-crack  interaction  region  for  a  nominal  strain 

rate  loading  of  10,000  s'1. 
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CHAPTER  7:  The  Effects  of  Grain  Boundary  Misorientations  on  Fracture  Propagation 
and  Crack  Stability  at  Crystalline-Amorphous  Energetic  Interfaces 

7.1  Introduction 

Local  failure  in  the  heterogeneous  crystalline-amorphous  structure  of  energetic 
materials,  the  defects  that  span  several  spatial  scales,  and  the  dynamic  thermo-mechanical 
loading  conditions  emphasize  that  a  fundamental  understanding  of  how  local  failure  occurs, 
at  the  interfaces  in  energetic  materials,  occurs  is  needed.  Therefore,  the  objective  of  this 
investigation  is  to  formulate  a  predictive  framework  that  accounts  for  the  interrelated 
microstructural  effects  of  crystalline  plasticity,  polymer  binder  relaxation,  crystalline- 
amorphous  fracture,  inherent  microstructural  defects  and  features,  and  dynamic  loading 
conditions  on  dynamic  crack  nucleation  and  growth  near  the  interfaces  in  energetic  materials. 
A  recently  developed  dislocation-density  based  crystalline  plasticity  [35]  approach,  which 
has  been  coupled  to  a  finite  viscoelasticity  formulation,  dynamic  fracture  formulation,  and 
dynamic  finite-element  techniques  will  be  used  to  understand  local  failure  at  RDX-estane 
interfaces. 

7.2  Results 

Computational  analyses  were  undertaken  to  investigate  and  predict  the  effects  of  low 
and  high  angle  GB  misorientations  on  crack  propagation  behavior  at  an  interface  for  a  RDX 
energetic  bi-crystal  with  a  viscoelastic  polymer  binder  between  the  crystals.  The  finite 
element  aggregate  consisted  of  2  grains  of  RDX  crystals  with  an  estane  binder  dispersed 
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between  the  crystals  for  an  overall  interface  model  of  1mm  x  1mm  (Figure  7.1)  with  a  pre¬ 
crack  length  of  50  pm  on  the  left  edge  of  the  model. 

Based  on  Gallagher  et  al.  [44],  the  RDX  grains  were  assumed  to  have  the  following 
three  independent  slip  systems  -  (010)  [001],  (021)  [100],  (02-1)  [100]  -  with  material 
properties  given  by  Annstrong  and  Elban  [1],  Annapragada  et  al.  [50], and  Clayton  and 
Becker  [60]  in  Table  7. 1 .  Clayton  and  Becker  [60]  give  a  range  of  yield  stress  for  RDX, 
from  approximately  160  to  640  MPa.  Therefore,  it  was  assumed  that  the  yield  stress  is  300 
MPa  for  this  study.  A  polymer  binder  of  estane  was  distributed  between  the  RDX  grains,  for 
a  spatially  uniform  interface.  Material  properties  for  the  estane  polymer  binder  were  obtained 
from  measurements  by  Mas  et  al.  [39]  and  Barua  and  Zhou  [23],  as  shown  in  Table  7. 1 . 

The  domain  of  influence  radius  for  fracture  nucleation  was  set  to  50  pm,  which  is  the 
width  of  the  estane  binder.  The  fracture  stress  was  set  to  five  times  the  yield  stress  of  RDX 
for  both  RDX  and  the  estane  binder.  Fixed  boundary  conditions  were  applied  to  the  model 
bottom  for  two-dimensional  plane-strain  loading  conditions.  A  tensile  velocity  was  applied 
on  the  top  surface,  as  shown  in  Figure  7.1,  for  dynamic  nominal  strain  rate  of  10,000  s’1,  with 
random  low  angle  GB  misorientations  of  less  than  20°  and  random  high  angle  GB 
misorientations  between  20°  and  70°. 

7.2.1  Bi-Crystal,  Random  Low  Angle  GB  Misorientations  for  a  Nominal  Strain  Rate  of 
10,000  s1  “ 

For  a  nominal  strain  rate  of  10,000  s’1,  at  a  nominal  strain  of  6.5%,  with  random  low 
angle  GB  misorientations  and  a  pre-existing  crack  of  length  50  pm,  the  normal  stress  (Figure 
7.2a)  was  a  maximum  at  the  interface  opposite  of  where  the  pre-crack  had  propagated  to  the 
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estane  binder  and  subsequently  tenninated.  The  maximum  nonnalized  normal  stress  was  4 
and  the  minimum  normalized  normal  stress  was  0  in  the  regions  that  had  unloaded  due  to 
crack  propagation  (stresses  were  normalized  by  the  static  yield  stress).  Another  crack 
nucleated  at  the  right,  bottom  free  boundary  of  the  model  due  to  high  stresses,  which  are 
most  likely  due  to  the  local  boundary  conditions.  The  lateral  stress  (Figure  7.2b)  had  a 
maximum  normalized  value  of  2,  which  occurred  where  the  secondary  crack  was  interacting 
with  the  estane  binder,  and  a  minimum  normalized  value  of  -0.5,  which  occurred  in  the 
binder  where  extensive  polymer  shearing  occurred  due  to  unloading  associated  with  crack 
propagation.  The  shear  stress  (Figure  7.2c)  had  a  maximum  nonnalized  value  of  1,  which 
occurred  in  the  region  of  crack  interaction  with  the  polymer  binder  for  the  secondary  crack, 
with  a  minimum  normalized  value  of  -0.2  occurring  in  the  region  dominated  by  shearing  due 
to  unloading  associated  with  crack  propagation.  The  pressure  (Figure  7. 2d)  was  a  maximum 
in  where  the  secondary  crack  interacted  with  the  estane  binder. 

To  further  understand  the  mechanics  of  crack  arrest  at  the  interface,  Figure  7.3  shows 
the  spatial  distribution  of  normal  stress,  shear  stress,  and  lateral  stress  for  the  thick,  dark  line 
shown  in  Figure  7. 2d.  The  normal  stress,  lateral  stress,  and  shear  stress  on  the  left  side  of  the 
RDX-estane  interface  (inside  of  the  left  RDX  crystal)  was  low,  with  a  maximum  normalized 
value  of  0.65  occurring  for  nonnal  stress.  As  the  spatial  location  moves  into  the  binder,  the 
nonnal  stress,  lateral  stress,  and  shear  stress  begin  to  rise.  At  the  right  RDX-estane  interface 
(at  the  inner  edge  of  the  right  RDX  crystal),  the  nonnal  stress  increases  sharply,  to  a 
maximum  normalized  value  of  4.25,  and  then  decreases  as  the  spatial  location  moves  away 
from  the  interface.  The  lateral  stress  decreases  as  the  spatial  location  moves  away  from  the 
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right  RDX-estane  interface.  Whereas,  the  shear  stress  increases  slightly  to  a  maximum 
nonnalized  value  of  approximately  0.6,  and  then  begins  to  decrease.  Therefore,  stresses  in 
the  binder  region  of  the  interface  are  lower  than  the  fracture  stress  due  to  the  viscous  nature 
of  the  polymer  binder.  This  further  substantiates  that  crack  propagation  was  arrested  due  to 
low  stresses  in  the  polymer  binder,  which  is  consistent  with  experimental  results  [29-3 1], 

The  temperature  (Figure  7.4a),  which  includes  the  effects  of  thermal  decomposition, 
plastic  work  heating,  viscous  dissipation,  and  thermal  conduction,  had  a  maximum 
nonnalized  value  of  1 .08  (temperatures  were  nonnalized  by  the  reference  temperature, 
293K),  which  occurred  in  the  shear  band  in  the  right  RDX  crystal.  The  shear  slip  (Figure 
7.4b),  which  is  the  plastic  slip  summed  over  all  slip  systems,  had  a  maximum  value  of  16%, 
while  small  amounts  of  shear  strain  localization  occuned  directly  on  the  opposing  side  of  the 
interface  from  the  crack  that  was  impinging  on  the  estane  binder.  The  dissipated  energy 
(Figure  7.4c),  which  includes  the  energy  dissipated  by  plasticity  and  viscous  dissipation,  had 
a  maximum  of  2E8  joules.  There  was  considerable  build  up  of  viscously  dissipated  energy  in 
the  binder  directly  ahead  of  the  crack  tip  as  a  result  of  the  impinging  crack-tip. 

The  mobile  dislocation  density  for  slip  system  one  (Figure  7.5a),  (010)  [001],  which 
was  the  most  active  mobile  slip  system,  had  a  maximum  nonnalized  value  of  450  (mobile 
dislocation  densities  were  nonnalized  by  the  initial  mobile  dislocation  density).  The 
immobile  dislocation  density  for  slip  system  one  (Figure  7.5b),  (010)  [001],  which  was  the 
most  active  immobile  slip  system,  had  a  maximum  nonnalized  value  of  45.  The  defonnation, 
for  random  low  angle  GB  misorientations,  was  dominated  by  slip  system  one,  (010)  [001]. 
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The  crack  tip  normalized  velocity  (velocity  was  normalized  by  the  dilatational  wave 
speed  of  3,379  m/s)  for  the  crack  that  propagated  from  the  pre-existing  crack  in  the  left  side 
of  the  model  is  shown  in  Figure  7.6  as  a  function  of  time.  The  crack  tip  normalized  velocity 
increases  slowly  until  a  time  of  approximately  3.1  psec,  where  the  nonnalized  velocity 
increased  sharply  to  a  peak  of  approximately  0.13  at  a  time  of  4.75  psec.  Peak  crack  tip 
velocity  was  attained  as  the  crack  approached  the  RDX-estane  interface,  where  the  velocity 
sharply  decreased  to  zero  at  the  interface  at  a  time  of  approximately  5.05  psec.  This  further 
substantiates  that  crack  arrest  occurred  as  the  crack  tip  approached  the  estane  binder. 

The  nominal  stress-strain  curve  (Figure  7.7)  increases  proportionally  with  strain  up  to 
a  nominal  strain  of  4%  and  a  nominal  normalized  stress  of  2.25  (stress  was  normalized  by  the 
static  yield  stress).  At  4%  nominal  strain,  the  stress  sharply  unloaded  briefly  due  to  crack 
propagation.  At  5%  nominal  strain,  the  stress  increased  with  strain  until  further  crack 
propagation  occurred  and  the  stress  unloaded.  The  after  cracks  propagated  through  both  sides 
of  the  bi-crystal  model,  at  a  nominal  strain  of  7%,  the  nominal  stress  oscillated  between  0  and 
1 .  Unloading  in  the  stress-strain  curve  was  affected  by  dynamic  crack  propagation  and 
nucleation  as  well  as  dynamic  plasticity. 

7.2.2  Bi-Crystal,  Random  High  Angle  GB  Misorientations  for  a  Nominal  Strain  Rate  of 
10,000  s1  ' 

For  a  nominal  strain  rate  of  10,000  s'1,  at  a  nominal  strain  of  7.5%  and  random  high 
angle  GB  misorientations,  the  nonnalized  nonnal  stress  (Figure  7.8a)  was  a  maximum  of  4 
and  the  minimum  was  0,  which  occuned  in  the  unloaded  regions  of  the  model,  due  to  crack 
nucleation  and  propagation.  Two  cracks  were  present  in  the  model;  one  propagated  from  the 
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pre-existing  crack  and  the  second  nucleated  at  the  bottom,  right  edge  of  the  model.  The  crack 
that  propagated  from  the  pre-existing  crack  had  two  sharp  regions,  which  are  associated  with 
a  change  in  the  cleavage  plane  that  had  maximum  stress.  One  sharp  change  in  crack  path 
occurred  approximately  150  pm  from  the  estane  binder,  while  the  second  occurred  at  the 
crack  tip  intersection  with  the  estane  binder.  This  is  a  result  of  the  high  random 
misorientations  that  affect  the  orientations  of  the  cleavage  normals,  and  further  substantiates 
the  importance  of  misorientations  in  aggregates,  as  compared  with  low  misorientations, 
where  smooth  crack  propagation  occurred.  The  lateral  stress  (Figure  7.8b)  had  a  maximum 
nonnalized  value  of  3  and  a  minimum  value  of  -1,  which  occurred  in  the  estane  binder, 
between  the  two  cracks.  The  shear  stress  (Figure  7.8c)  had  a  maximum  normalized  value  of 
1.5  and  a  minimum  value  of  -0.5  that  occurred  just  above  the  crack  tip  region.  The  pressure 
(Figure  7.8d)  had  a  maximum  normalized  value  of  3  and  a  minimum  value  of  -0.5. 

Maximum  temperature  accumulation  (Figure  7.9a)  occurred  at  the  two  sites  where  the 
cleavage  planes  changed  with  a  maximum  nonnalized  value  of  1.3.  The  shear  slip  (Figure 
7.9b)  was  a  maximum  in  the  first  sharp  change  in  crack  path,  with  a  maximum  value  of  90%. 
The  maximum  dissipated  energy  (Figure  7.9c)  was  4E08  joules  and  occurred  in  the  same 
region  as  the  maximum  temperature  accumulation. 

The  immobile  dislocation  density  buildup  for  slip  system  one  (Figure  7. 10a)  had  a 
maximum  nonnalized  value  of  40,  while  slip  system  two  (Figure  7. 10b)  had  a  maximum 
nonnalized  accumulation  of  600.  Slip  system  three  (Figure  7.10c)  had  a  maximum 
nonnalized  buildup  of  50.  The  maximum  buildup  of  dislocation  densities  and  buildup  for  the 
low  angle  misorientations  was  an  order  of  magnitude  lower  than  the  high  angle 
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misorientations.  This  demonstrates  that  high  angle  misorientations  facilitate  more  inelastic 
deformation  for  RDX  crystals  than  low  angle  misorientations.  This  is  consistent  with 
experimental  observations  that  defonnation  is  highly  anisotropic  and  orientation  dependent 
for  RDX  crystals  [57,  61]  due  to  the  low  number  of  unique  slip  systems  available  for 
inelastic  deformation. 

To  understand  the  sharp  pointed  crack  path  in  the  left  RDX  crystal,  Figure  7.1 1  shows 
the  temporal  evolution  of  normal  stress,  lateral  stress,  pressure,  and  temperature  for  the 
critical  element  region  shown  in  Figure  7.10c  with  the  thick,  dark  circle.  Nonnal  stress, 
pressure,  and  lateral  stress  all  increase  up  to  a  time  of  4.5  psec,  at  which  time,  they  sharply 
increased  to  a  nonnalized  value  of  10,  6,  and  4.75,  respectively,  which  was  a  direct 
consequence  of  the  resulting  crack  propagating  near  the  critical  element  region  and  the 
subsequent  changing  the  cleavage  plane  of  crack  propagation  direction  due  to  lattice  rotations 
as  a  result  of  plastic  deformation.  The  temperature  remained  at  the  initial  temperature  until 
approximately  4.5  psec,  where  the  temperature  sharply  increased  to  a  normalized  value  of 
1 .3 1  at  a  time  of  7  psec.  The  temperature  then  gradually  decreased  to  the  effects  of  thermal 
conduction  for  the  remaining  time.  Therefore,  lattice  rotation,  due  to  deformation  of  high 
angle  GB  misorientations,  raises  the  probability  of  switching  cleavage  plane  propagation 
directions,  which  can  result  in  large  temperature  buildups  due  to  the  temporary  pinning  of  the 
crack  tip. 

The  crack  tip  velocity  for  the  crack  in  the  left  RDX  crystal  is  shown  in  Figure  7.12. 
The  crack  tip  velocity  slowly  increases  until  a  time  of  3  psec,  at  which  time,  the  crack  tip 
nonnalized  velocity  peaked  at  approximately  0.08  at  a  time  of  5  psec.  This  peak  in  crack 
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velocity  corresponds  to  the  sharp  increase  in  nonnal  stress,  lateral  stress,  pressure,  and 
temperature  shown  in  Figure  7.1 1  which  is  associated  with  the  subsequent  arresting  of  the 
crack  which  is  evident  due  to  the  sudden  drop  in  crack  tip  velocity  after  the  first  peak  is 
attained.  The  crack  tip  velocity  remains  low  until  a  time  of  approximately  6.75  psec,  at  which 
time  the  crack  tip  velocity  becomes  unbounded  which  is  associated  with  unstable  crack 
growth.  The  crack  tip  propagated  to  the  interface  where  it  was  arrested  due  to  low  stresses  in 
the  binder  regions.  Therefore,  the  crack  propagated  to  the  critical  element  region  shown  in 
Figure  7. 10c,  and  was  arrested  for  approximately  2  nsec  due  to  a  change  in  cleavage  plane 
propagation  direction.  At  which  point  the  stress  on  the  new  cleavage  plane  built  up.  At  a  time 
of  6.75  psec,  the  stresses  on  the  new  cleavage  plane  violated  the  failure  criterion  and 
subsequently  propagated  unstably  to  the  RDX-estane  interface  where  it  was  arrested.  This  is 
in  comparison  to  the  low  angle  GB  misorientations  case,  where  the  crack  velocity 
continuously  built  up  until  the  crack  tip  reached  the  RDX-estane  interface,  where  it  was 
arrested. 

The  nominal  stress-strain  behavior  for  low  and  high  angle  misorientations  is  shown  in 
Figure  7.13.  For  high  angle  misorientations,  the  stress  increased  up  to  a  nominal  strain  of  4%, 
and  a  nominal  stress  of  2.4.  At  4%  nominal  strain,  the  pre-existing  crack  propagated,  which 
is  marked  by  sharp  drop  in  stress  at  5%  nominal  strain.  The  stress  began  to  increase  up  to  a 
strain  of  7%,  which  indicates  the  crack  has  been  arrested.  Then,  the  crack  continued  to 
propagate,  which  is  marked  by  a  sharp  drop  in  stress  at  a  nominal  strain  of  7.5%.  The  stress- 
strain  behavior  (Figure  7.13)  demonstrates  that  high  angle  misorientations  lead  to  a  larger 
build  up  of  stress  due  to  crack  arrest  as  well  as  large  regions  of  plasticity  in  comparison  to  the 
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low  angle  misorientation  case.  For  low  angle  misorientations,  the  bi-crystal  unloads  due  to 
fracture  at  an  earlier  strain  in  comparison  to  high  angle  misorientations.  The  peak  stress  for 
high  angle  misorientations  is  higher  than  that  of  low  angle  misorientations,  which  is  due  to 
the  temporary  arresting  of  the  crack. 

7.2.3  16  Crystal  Aggregate,  Random  Low  Angle  GB  Misorientations  for  a  Nominal  Strain 
Rate  of  10,000  s' 1 

The  finite-element  aggregate  of  16  grains  of  RDX  crystals  with  an  average  grain  size 
of  250  pm  for  an  overall  aggregate  size  of  1mm  x  1mm  (Figure  7.14)  consisted  of 
approximately  7,000  quadrilateral  elements.  The  RDX  crystalline  morphology  was  generated 
utilizing  a  Voronoi  tessellation  scheme.  Based  on  Gallagher  et  al.  [44],  the  RDX  grains  were 
assumed  to  have  the  following  three  independent  slip  systems  -  (010)  [001],  (021)  [100], 
(02-1)  [100]  -  with  material  properties  given  by  Armstrong  and  Elban  [1]  and  Annapragada 
et  al.  [50]  in  Table  7. 1 .  A  polymer  binder  of  estane,  with  a  volume  fraction  of  approximately 
15%,  was  randomly  distributed  between  the  RDX  grains.  The  material  properties  for  the 
estane  polymer  binder  were  based  on  measurements  by  Mas  et  al.  [39]  and  Barua  and  Zhou 
[23],  as  shown  in  Table  7. 1 .  A  pre-existing  crack  of  100  pm  was  used  in  the  left  side  of  the 
model. 

The  pre-existing  crack  was  placed  to  model  a  crack  that  has  already  nucleated.  The 
immobile  dislocation  density  for  slip  system  one  (Figure  7.15a),  (010)  [001],  for  a  nominal 
strain  rate  of  10,000  s'1,  at  a  nominal  strain  of  10%,  accumulated  a  maximum  dislocation 
density  of  180  (dislocation  densities  have  been  normalized  by  the  initial  densities).  The 
dislocation  density,  which  is  a  measure  of  inelastic  activity  on  slip  system  one,  had  localized 
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at  the  crack  tip  and  the  crystal  edges  where  shear  deformation  occurred  in  the  polymer 
binder.  Shear  slip  (Figure  7.15b)  built  up  at  the  crack  tip  and  the  crystal  facets  associated 
with  shear  defonnation  in  the  polymer  binder.  The  energy  dissipated  due  to  viscous 
dissipation  and  plasticity  (Figure  7.15c)  localized  in  the  binder  regions,  where  viscous 
shearing  dominated  the  defonnation,  and  at  the  crack  tip  and  crystal  facets. 

The  pre-crack  in  the  aggregate,  did  not  propagate,  but  localized  plasticity  did  fonn  at 
the  tip  of  the  crack  tip.  Shear  deformation  in  the  polymer  allowed  relaxation  of  the  binder, 
which  constrained  the  crack  from  propagating,  due  to  limiting  stresses  to  below  the  fracture 
stress.  A  combination  of  polymer  shear  relaxation  and  plasticity  at  the  crack  tip  prevented  the 
stresses  on  the  cleavage  planes  from  violating  the  failure  criterion.  Therefore,  polymer  binder 
relaxation  and  crack  tip  plasticity  arrested  the  crack.  These  results,  for  an  aggregate  with 
crystal  morphology,  substantiate  that,  although,  at  the  interfacial  scale,  cracks  can  propagate 
and  blunt  due  to  plasticity,  but  also  cracks  can  be  arrested  due  to  the  relaxation  associated 
with  the  polymer  binder. 

7.3  Summary 

The  effects,  of  GB  misorientations  on  crack  propagation  towards  a  polymer  binder  in 
a  bi-crystalline  model,  have  been  studied  at  a  nominal  strain  rate  of  10,000  s'1.  Pre-existing 
cracks  propagated  towards  the  binder  and  were  arrested  due  to  low  stresses  in  the  binder.  The 
low  binder  stresses  are  a  direct  consequence  of  the  polymer  material,  which  illuminates  why 
cracks  propagate  in  crystalline  regions  and  tenninate  at  the  binder. 

For  low  angle  GB  misorientations,  the  pre-existing  crack  propagated  in  a  straight  line 
towards  the  estane  binder.  The  crack  tip  velocity  increased  continuously  until  the  crack 
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reached  the  estane  binder,  where  it  was  arrested.  The  inelastic  defonnation  was  dominated  by 
one  slip  system.  For  high  angle  misorientations,  the  pre-existing  crack  propagated  towards 
the  estane  binder,  but  was  arrested  briefly  while  the  crack  propagation  direction  changed  from 
one  cleavage  plane  to  another  due  to  the  high  misorientations  and  further  exacerbated  by 
lattice  rotations  due  to  inelastic  defonnation.  During  this  brief  period  of  crack  arrest,  the 
inelastic  deformation  built  up  on  three  slip  systems  and  the  temperature  increased  at  the  crack 
tip.  The  stress  on  the  new  cleavage  plane  attained  the  failure  stress  and  the  crack  unstably 
propagated  towards  the  binder,  where  it  was  arrested.  For  an  aggregate,  the  relaxation 
associated  with  the  polymer  binder  can  result  in  crack  arrest  due  to  the  lowering  of  stresses  on 
the  cleavage  planes  by  relaxing  mechanical  strains  through  viscous  dissipation  in  the 
polymer. 

Results  from  this  study  indicate  that  the  interrelated  effects  of  dynamic  strain  rate 
loading,  polymer  binder  response,  GB  misorientation,  and  pre-existing  cracks  significantly 
affect  crack  propagation  towards  an  interface.  Pre-existing  cracks  propagate  differently 
depending  on  the  level  of  misorientation.  Plasticity  in  RDX  is  dominated  by  different  slip 
systems  due  to  the  highly  anisotropic  orientation  dependent,  inelastic  deformation  behavior. 
Cracks  can  be  arrested  due  to  plasticity  and  polymer  binder  relaxation. 
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7.4  Tables  and  Figures 


Table  7.1:  Material  Properties  for  RDX  and  Estane 


RDX 

Estane  Binder 

Density  (!  ) 

1 820 

TT90 

I  ! 

Yield  Stress  (MPa) 

300 

— 

Elastic  Modulus  (GPa) 

18.4 

— 

Poisson  Ratio 

0.22 

.499 

! 

Thermal  Conductivity  (  '  ) 

0.29 

0.14 

f 

Specific  Heat  (cp) 

!  "  ! 

1260 

1500 

Rate  Sensitivity  Parameter 

.01 

— 

Fracture  Stress  (GPa) 

1.5 

1.5 

Initial  Mobile  Dislocation 

Density  (m"  ) 

lelO 

— 

Initial  Immobile  Dislocation 
Density  (m"2) 

lel2 

— 
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Figure  7. 1 :  The  RDX-estane  bi-crystal  with  loading,  boundary  conditions,  and  grain 

morphology. 
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Figure  7.2:  (a)  Normal  Stress,  (b)  lateral  stress,  (c)  shear  stress,  and  (d)  pressure  for  a 
nominal  strain  rate  loading  of  10,000  s'1,  at  a  nominal  strain  of  6.5%,  and  random  low  angle 
GB  misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  7.3:  Normal  stress,  shear  stress,  and  lateral  stress  spatial  distribution  across  the  RDX- 
estane  interface  for  a  nominal  strain  rate  loading  of  10,000  s'1,  at  a  nominal  strain  of  6.5%, 

and  random  low  angle  GB  misorientations. 
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Figure  7.4:  (a)  Temperature,  (b)  shear  slip,  and  (c)  dissipated  energy  for  a  nominal  strain  rate 
loading  of  10,000  s’1,  at  a  nominal  strain  of  6.5%,  and  random  low  angle  GB  misorientations. 

(Spatial  dimensions  are  in  units  of  meters) 
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Figure  7.5:  (a)  Mobile  dislocation  density  and  (b)  immobile  dislocation  density  for  slip 
system  one,  which  was  the  most  active  slip  system,  for  a  nominal  strain  rate  loading  of 
10,000  s’1,  at  a  nominal  strain  of  6.5%,  and  random  low  angle  GB  misorientations.  (Spatial 

dimensions  are  in  units  of  meters) 
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Velocity  (normalized  by  dilatational  wave  speed) 


Figure  7.6:  Crack  tip  velocity  versus  time  for  the  crack  propagating  from  the  initial  pre-crack, 
for  a  nominal  strain  rate  loading  of  10,000  s'1  and  random  low  angle  GB  misorientations. 
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Figure  7.7:  Nominal  stress-strain  curve  for  a  nominal  strain  rate  loading  of  10,000  s'1  and 

random  low  angle  GB  misorientations. 
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Figure  7.8:  (a)  Normal  Stress,  (b)  lateral  stress,  (c)  shear  stress,  and  (d)  pressure  for  a 
nominal  strain  rate  loading  of  10,000  s'1,  at  a  nominal  strain  of  7.5%,  and  random  high  angle 
GB  misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  7.9:  (a)  Temperature,  (b)  shear  slip,  and  (c)  dissipated  energy  for  a  nominal  strain  rate 
loading  of  10,000  s'1,  at  a  nominal  strain  of  7.5%,  and  random  high  angle  GB 
misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  7.10:  Immobile  dislocation  densities  for  (a)  slip  system  one,  (b)  slip  system  two,  and 
(c)  slip  system  three  for  a  nominal  strain  rate  loading  of  10,000  s  ,  at  a  nominal  strain  of 
7.5%,  and  random  high  angle  GB  misorientations.  (Spatial  dimensions  are  in  units  of  meters) 
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Figure  7.11:  Time  history  plot  of  normal  stress,  pressure,  lateral  stress,  and  temperature  for  a 
nominal  strain  rate  loading  of  10,000  s'1  and  random  high  angle  GB  misorientations  in  the 

critical  element  location. 
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Figure  7.12:  Crack  tip  velocity  versus  time  for  the  crack  propagating  from  the  initial  pre¬ 
crack,  for  a  nominal  strain  rate  loading  of  10,000  s'1  and  random  high  angle  GB 

misorientations. 
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Figure  7.13:  Stress  versus  strain  curves  for  a  nominal  strain  rate  loading  of  10,000  s  ,  for 
random  low  and  high  angle  GB  misorientations. 
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Figure  7. 14:  The  RDX-estane  aggregate  with  loading  and  boundary  conditions,  pre-existing 

crack,  grain  morphology. 
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Figure  7.15:  (a)  Immobile  dislocation  densities  for  slip  system  one,  (b)  inelastic  shear  slip, 
and  (c)  dissipated  energy  due  to  viscous  effects  and  plasticity  for  a  nominal  strain  rate 
loading  of  10,000  s'1,  at  a  nominal  strain  of  10%,  and  random  low  angle  GB  misorientations. 

(Spatial  dimensions  are  in  units  of  meters) 
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CONCLUSIONS 


Through  these  analyses,  high  dynamic  pressure  loads  and  strain  rate  loading 
conditions  have  been  investigated  to  understand  the  effects  of  unbounded  local 
decomposition  temperatures  and  hot  spot  formation  at  the  microscale.  This  work  showed  that 
hot  spot  fonnation  is  coupled  to  the  inherent  crystal  defects  and  how  defects  result  in  energy 
localization  at  defect  sites.  This  work  also  showed  that  the  microstructural  response  of 
energetic  materials  is  a  direct  function  of  the  highly  anisotropic  inelastic  deformation  that  is 
unique  to  the  energetic  crystal  RDX,  which  is  a  direct  consequence  of  the  limited  slip 
systems  for  inelastic  deformation. 

For  a  polymer,  where  the  operating  temperature  was  below  the  glass  transition 
temperature,  hot  spot  formation  resulted  at  the  peripheries  of  the  voids  where  inelastic 
deformation  was  dominant.  However,  for  a  polymer  where  the  operating  temperature  was 
above  the  glass  transition  temperature,  hot  spot  formation  occurred  at  the  peripheries  of  the 
crystals  where  shearing  in  the  polymer  region  was  dominant  and  resulted  in  crystal  to  crystal 
interaction.  Crystal  to  crystal  interactions  occurred  due  a  thin  layer  of  polymer  that  was 
hydro-statically  constrained  between  the  crystals  that  facilitated  crystal  sliding  and  allowed 
for  the  transfer  of  loads  from  one  crystal  to  another.  This  transfer  of  load,  when  large, 
resulted  in  inelastic  deformation  at  the  crystal  peripheries  and  subsequently  hot  spot 
formation. 

This  work  also  investigated  the  effects  of  fracture  at  the  microscale,  where  cracks 
nucleated  from  internal  void  porosity  due  to  large  stresses.  Crack  nucleation  on  the  cleavage 
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planes  at  the  peripheries  of  the  voids  was  subsequently  followed  by  crack  propagation 
towards  the  binder  where  crack  arrest  occurred.  Crack  arrest  at  the  binder  occurred  due  to  the 
viscous  nature  of  the  polymer  material,  which  resulted  in  low  stress  buildup  in  the  polymer 
binder  regions,  which  acted  to  arrest  the  crack.  Where  cracks  impinged  upon  the  interface  of 
the  RDX-binder,  a  buildup  of  inelastic  defonnation  and  subsequent  temperature  accumulation 
occurred,  which  resulted  in  hot  spot  formation  at  the  interface  due  to  the  effects  of  fracture. 

These  analyses  have  identified  the  dominant  microstructural  mechanisms  that  affect 
hot  spot  fonnation  and  failure  evolution  in  energetic  materials,  and  the  results  from  these 
analyses  can  be  used  to  tailor  energetic  material  aggregates  to  a  desired  response  based  on 
the  microstructural  constituents,  defects,  and  behavior. 
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FUTURE  RESEARCH  RECOMMENDATIONS 


1)  Develop  a  crack  nucleation  model  based  on  dislocation  density  pileup  release 
relaxation  at  the  RDX-binder  interface,  to  allow  for  the  relaxation  of  dislocation 
density  pileups,  which  can  nucleate  micro-cracks  at  the  grain  boundaries. 

2)  Characterize  void-to-void  interactions  to  account  for  void  coalescence  and  shear 
strain  localization  in  the  intervoid  regions. 

3)  Investigate  the  effects  of  debonding  on  hot  spot  formation  and  local  and  global  failure 
of  energetic  aggregates. 

4)  Investigate  shear  fracture  under  compressive  loadings,  and  crystal-to-crystal 
interactions  under  compression  utilizing  a  contact  algorithm. 

5)  Link  thermal  decomposition  to  chemical  and  combustion  reactions,  to  use  atomistic 
and  molecular  dynamic  models  to  accurately  predict  the  decomposition  behavior  of 
energetic  materials. 

6)  Use  these  predictions  to  improve  equations  of  state  for  shock  loadings  of  energetic 
aggregate  materials. 
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